{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fault detection\n",
    "\n",
    "We'll consider a problem of identifying faults that have occurred in a system based on sensor measurements of system performance.\n",
    "\n",
    "# Topic references\n",
    "- [Samar, Sikandar, Dimitry Gorinevsky, and Stephen Boyd. \"Likelihood Bounds for Constrained Estimation with Uncertainty.\" Decision and Control, 2005 and 2005 European Control Conference. CDC-ECC'05. 44th IEEE Conference on. IEEE, 2005.](http://web.stanford.edu/~boyd/papers/pdf/map_bounds.pdf)\n",
    "\n",
    "# Problem statement\n",
    "\n",
    "Each of $n$ possible faults occurs independently with probability $p$.\n",
    "The vector $x \\in \\lbrace 0,1 \\rbrace^{n}$ encodes the fault occurrences, with $x_i = 1$ indicating that fault $i$ has occurred.\n",
    "System performance is measured by $m$ sensors.\n",
    "The sensor output is\n",
    "\\begin{equation}\n",
    "y = Ax + v = \\sum_{i=1}^n a_i x_i + v,\n",
    "\\end{equation}\n",
    "where $A \\in \\mathbf{R}^{m \\times n}$ is the sensing matrix with column $a_i$ being the **fault signature** of fault $i$,\n",
    "and $v \\in \\mathbf{R}^m$ is a noise vector where $v_j$ is Gaussian with mean 0 and variance $\\sigma^2$.\n",
    "\n",
    "The objective is to guess $x$ (which faults have occurred) given $y$ (sensor measurements).\n",
    "\n",
    "We are interested in the setting where $n > m$, that is, when we have more possible faults than measurements.\n",
    "In this setting, we can expect a good recovery when the vector $x$ is sparse.\n",
    "This is the subject of compressed sensing.\n",
    "\n",
    "# Solution approach\n",
    "To identify the faults, one reasonable approach is to choose $x \\in \\lbrace 0,1 \\rbrace^{n}$ to minimize the negative log-likelihood function\n",
    "\n",
    "\\begin{equation}\n",
    "\\ell(x) = \\frac{1}{2 \\sigma^2} \\|Ax-y\\|_2^2 +  \\log(1/p-1)\\mathbf{1}^T x + c.\n",
    "\\end{equation}\n",
    "\n",
    "However, this problem is nonconvex and NP-hard, due to the constraint that $x$ must be Boolean.\n",
    "\n",
    "To make this problem tractable, we can relax the Boolean constraints and instead constrain $x_i \\in [0,1]$.\n",
    "\n",
    "The optimization problem\n",
    "\n",
    "\\begin{array}{ll}\n",
    "\\mbox{minimize} &  \\|Ax-y\\|_2^2 + 2 \\sigma^2 \\log(1/p-1)\\mathbf{1}^T x\\\\\n",
    "\\mbox{subject to} &  0 \\leq x_i \\leq 1, \\quad i=1, \\ldots n\n",
    "\\end{array}\n",
    "\n",
    "is convex.\n",
    "We'll refer to the solution of the convex problem as the **relaxed ML** estimate.\n",
    "\n",
    "By taking the relaxed ML estimate of $x$ and rounding the entries to the nearest of 0 or 1, we recover a Boolean estimate of the fault occurrences.\n",
    "\n",
    "# Example\n",
    "\n",
    "We'll generate an example with $n = 2000$ possible faults, $m = 200$ measurements, and fault probability $p = 0.01$.\n",
    "We'll choose $\\sigma^2$ so that the signal-to-noise ratio is 5.\n",
    "That is,\n",
    "\\begin{equation}\n",
    "\\sqrt{\\frac{\\mathbf{E}\\|Ax \\|^2_2}{\\mathbf{E} \\|v\\|_2^2}} = 5.\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "np.random.seed(1)\n",
    "\n",
    "n = 2000\n",
    "m = 200\n",
    "p = 0.01\n",
    "snr = 5\n",
    "\n",
    "sigma = np.sqrt(p*n/(snr**2))\n",
    "A = np.random.randn(m,n)\n",
    "\n",
    "x_true = (np.random.rand(n) <= p).astype(np.int)\n",
    "v = sigma*np.random.randn(m)\n",
    "\n",
    "y = A.dot(x_true) + v"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below, we show $x$, $Ax$ and the noise $v$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x11ae42518>]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAGMNJREFUeJzt3X2QXXV9x/HPl4TwFJ6zQUyCCbq2xrYWuoM4VGsLQsLUpD42mbZSS00djQ8DVnHoUIf+0RGnxTpGMFZGZXgKqDU6i+hQLIiA2ZAHCCGwBDDL5mHJI5Bsspt8+8c9u3tz9z6ce/fc8/Dj/ZrZyb3nnnvPd3/33E/OOd9z9pq7CwAQlmOyLgAAkDzCHQACRLgDQIAIdwAIEOEOAAEi3AEgQIQ7AASIcAeAABHuABCgyVkteNq0aT579uysFg8AhbR69eqX3b2j0XyZhfvs2bPV09OT1eIBoJDM7MU483FYBgACRLgDQIAIdwAIEOEOAAEi3AEgQA3D3cxuMbMdZvZkjcfNzL5hZr1mtt7Mzk++TABAM+JsuX9P0rw6j8+X1Bn9LJF008TLAgBMRMNwd/cHJe2qM8tCST/wkkclnWZmZydV4ESsfnG3Nm7dF3v+n63v1579h1pa1r7BIf1k7Us1H3d33bO6T4NDh0endT+xVbtei7e8waHD+uHqPk30axFXvbBLm7a9UvWxrXsP6P6N22s+1931o8f7tP/QsJ7q36fHf7dbkrT6xV16elv8cW5k7/4h/XRd/1HTHtu8U89ur153s57Z/op++/z4VXr48BGtWLVFR47k46snV67r194DQ7HmPXLEtaJni4YOH2lzVdlb37dH6/v2VH1s12uH1P3E1rbXsPfA+HU0b5I45j5D0pay+33RtHHMbImZ9ZhZz8DAQAKLru9DN/1G8//roVjz9u3er6W3r9HS29e0tKwv3bNen7tzbc2Q+9WmAX3h7nW64eebJEk7XhnUp257XJ+8dXWs1//37o26+u51+nXvyy3VN+IjNz+iy77+YNXHPrDsN7ry+7UvLPvt87t01Yp1uv6nT+nybzykD37rN5KkD930iOZ9Pd44x/H5u9boM3es0fMvvzY67a+XP6r33Vi97mZdeuOD+ui3Hxk3/ZaHn9cXf7hed/VsqfKsdD038Ko+e8caXb1ibaz5f7TmJX3xnvVa/uDmNleWvQXffFgLvvlw1cf+6dYefeq2x7XjlcG21nDVXWv1mTvWaPPAq21dzkQkEe5WZVrVTR93X+7uXe7e1dHR8OrZVB0cLm3x9O850NLz+/eWVqYDhw5XfXzfYGkLbODVg5KkQ9HyXoq5vO37Ss97dXC4pfri2Lav/gfitUPDUS3t/eD07ym9/sHh6mPZLjujvag9++NtLbfTyHo0MhaNjOxx7ny1tT3PUPTtLn2ehg+3d+9r5HM7OJTfPaUkwr1P0qyy+zMl5Xt/BQACl0S4r5T0seismQsl7XX39h/0AgDU1PAPh5nZHZLeK2mamfVJ+ldJx0qSu98sqVvS5ZJ6Je2X9PF2FQsAiKdhuLv74gaPu6RPJ1YRAGDCuEIVAAJEuANAgAh3AAgQ4Q4AASLcASBAhDsABIhwB4AAEe4AECDCHQACRLgDQIAIdwAIEOEOAAEi3AEgQIQ7AASIcAeAABHuABAgwh0AAkS4A0CACHcACBDhDgABItwxTuk7zwEUGeEOAAEi3AEgQIQ7AASIcAeAABHuGId+KlB8hDsABIhwB4AAEe4AECDCHQACFCvczWyemW0ys14zu6bK4+eY2QNmtsbM1pvZ5cmXirTQTwWKr2G4m9kkScskzZc0V9JiM5tbMdu/SFrh7udJWiTpW0kXCgCIL86W+wWSet19s7sfknSnpIUV87ikU6Lbp0rqT65EAECzJseYZ4akLWX3+yS9s2Ker0j6hZl9RtJJki5JpDoAQEvibLlblWmVh2UXS/qeu8+UdLmkW81s3Gub2RIz6zGznoGBgearBQDEEifc+yTNKrs/U+MPu1wpaYUkufsjko6XNK3yhdx9ubt3uXtXR0dHaxWj7fiTv0DxxQn3VZI6zWyOmU1RqWG6smKe30m6WJLM7G0qhTub5gCQkYbh7u7DkpZKuk/SRpXOitlgZteb2YJotqslfcLM1km6Q9LfO5t/AJCZOA1VuXu3pO6KadeV3X5K0kXJlgYAaBVXqAJAgAh3jMPxNKD4CHcACBDhDgABItwBIECEOwAEiHBPmFm1v9ZQezpqs6p/+QLVsH6hEuGesFrXbhXpmq68lOqctxNbkdYvpINwB4AAEe4AECDCHQACRLgnjIZqcmioxsf6hUqEe8KCaKjmpJGZlzqKoEjrF9JBuANAgAh3AAgQ4R4Z2atteec25m7xyO7z6PLiPi+qLMud7wmPUdzlZPVbjv5++TnE0Wwleao9C2mto6PLy/F4E+4AECDCPTJyskHL5xzEPFth5KyG0eXFfV5UWRrnRNTamZjwGMWU2Vkyo79ffs48abaSPNWehbTW0dHl5Xi8CXcACBDhDgABItwjNFRj1EBDNXU0VJtDQ3UM4Q4AASLcIzRUY9RAQzV1NFSbQ0N1DOEOAAEi3AEgQIR7hIZqjBpoqKaOhmpzaKiOIdwBIECEe4SG6hiuUM1Pk4yGanNoqI4h3AEgQIQ7AASIcI/QUI1RAw3V1NFQbQ4N1TGxwt3M5pnZJjPrNbNraszzUTN7ysw2mNntyZYJAGjG5EYzmNkkScskvU9Sn6RVZrbS3Z8qm6dT0pclXeTuu81sersKbhcaqmNqbY3QUE0fDdXmpP094Xke7zhb7hdI6nX3ze5+SNKdkhZWzPMJScvcfbckufuOZMsEADQjTrjPkLSl7H5fNK3cWyW91cweNrNHzWxetRcysyVm1mNmPQMDA61VDABoKE64V9vvqNxvnyypU9J7JS2W9N9mdtq4J7kvd/cud+/q6OhotlYAQExxwr1P0qyy+zMl9VeZ5yfuPuTuz0vapFLYo4aYJ8kAaEHcs9BCFifcV0nqNLM5ZjZF0iJJKyvm+R9Jfy5JZjZNpcM0m5MstChqNUjjNk4xJs/Nqrxh/UKlhuHu7sOSlkq6T9JGSSvcfYOZXW9mC6LZ7pO008yekvSApH92953tKjrPam0xFGlLIi+l5vkc4rwp0vqFdDQ8FVKS3L1bUnfFtOvKbrukq6IfAEDGuEIVAAJEuAMIDkepCPfE0VBNDg3V+Fi/UIlwT1gQDdWsC4jQUI2vSOsX0kG4A0CACHcACBDhnhH2ogG0E+GeMBqqyaGhGh/rFyoR7gkLoqGak1ppqMaXl/cM+UG4A0CACHcACBDhnhEOOQDtw1Eqwj1xNFSTQ0M1PtYvVCLcExZEQzXrAiLs3cRXpPUL6SDcASBAhDsABIhwj6S9V8teNNA+HNIj3BNHQzU5NFTjY/1CJcI9YUE0VHNSKltf8RVp/UI6CHcACBDhDgABItxHpbtby0400D4cpSLcE0dDNTk0VONj/UIlwh0AAkS4JyyEs2XycsyIs2XiK9T6hVQQ7gAQIMI9kv4VqmxpAe3Cp4twTxwN1eTQUI2P9QuVCHcACBDhnrAQGqp5aWTmpY4iKNL6hXTECnczm2dmm8ys18yuqTPfh83MzawruRIBAM1qGO5mNknSMknzJc2VtNjM5laZ72RJn5X0WNJFpiHt7R62s4D2YU8m3pb7BZJ63X2zux+SdKekhVXm+zdJN0gaTLC+wqGhiiywfqFSnHCfIWlL2f2+aNooMztP0ix3/1mCtQEAWhQn3KttEozu85jZMZJulHR1wxcyW2JmPWbWMzAwEL/KAgmioZqTUvNSRxEUaf1COuKEe5+kWWX3Z0rqL7t/sqQ/kPQrM3tB0oWSVlZrqrr7cnfvcveujo6O1qsGANQVJ9xXSeo0szlmNkXSIkkrRx50973uPs3dZ7v7bEmPSlrg7j1tqbhN+A5VIBxpfbzyfLpuw3B392FJSyXdJ2mjpBXuvsHMrjezBe0usGhoqCILrF+oNDnOTO7eLam7Ytp1NeZ978TLAgBMBFeoJiyIhmrWBUQKNGSZK9L6hXQQ7gAQIMI9kn5jhC0toF3S2pHJ8w4T4Z4wGqrIAusXKhHuABAgwj1hQTRUc1Jrns8hzpu8vGfID8IdAAJEuEdGNnxa3v6JueU0soU1ury4z4sqy3L7bMJjFHc5Wf2Wo79f9lvBrY51HmrP0tjH6fU9DhLhDgBBItwrtHzOQXS2QqOzZSofj3uWw8iXRWd5TsRIqe2uIbMvxh79/fJz5kncSkbXrxzVngVOGhpDuCcsiIZq1gVECjRkmSvS+oV0EO4AECDCPZJ6Q7XifsPn0VBtvzw1VFt8v/NQe5bGTlRId3l5RLgDQIAI9wo0VOvUQEM1dTRUm0NDdQzhDgABItwTFsTZMsUpFZEirV9IB+EemXDDkitUk1sODVWuUG1RWuvo6PJyPN6EOwAEiHCvMNGGau2HaajGXg4N1VE0VJtDQ3UM4Q4AASLcE1brCFyRGl55OY5YoCHLXJHWL6SDcI+kf4WqH3W/4fNoqLZfnhqqFf/Gf172tWeJK1THEO4AECDCvQIN1To10FBNHQ3V5tBQHUO4A0CACPeE1ToGV6iGV05Kfb0fP25GodYvpIJwr8AVqnVqoKGamso/DR37eTmoPUue8nuY59Em3AEgQIR7BRqqdWqgoZo6GqrNoaE6hnAHgADFCnczm2dmm8ys18yuqfL4VWb2lJmtN7P7zexNyZdaFAH8yd+sC4gUaMgyV6T1C+loGO5mNknSMknzJc2VtNjM5lbMtkZSl7v/kaR7JN2QdKHtlvoVqjRUay+HhipXqLYo/StU8zvecbbcL5DU6+6b3f2QpDslLSyfwd0fcPf90d1HJc1MtkwAQDPihPsMSVvK7vdF02q5UtK91R4wsyVm1mNmPQMDA/GrTBEN1To10FBNHQ3V5tBQHRMn3KsNV9V9ETP7W0ldkr5W7XF3X+7uXe7e1dHREb9KAEBTJseYp0/SrLL7MyX1V85kZpdIulbSn7n7wWTKK54QrlDNS6k5KaMQirR+IR1xttxXSeo0szlmNkXSIkkry2cws/MkfVvSAnffkXyZ7Zf6d6jyJ39rL4eGKt+h2qLUG6rpLKYlDcPd3YclLZV0n6SNkla4+wYzu97MFkSzfU3SVEl3m9laM1tZ4+UAACmIc1hG7t4tqbti2nVlty9JuK7M0FCtUwMN1dTRUG0ODdUxXKEKAAEi3AEgQIR7JKkrVBt9QXYRrlCt1ZRLu6Ga+hkgOWqoqsn3u7JR/3qV+p/8zfFwE+4AECDCvQIN1To10FBNHQ3V5tBQHUO4A0CACHcACBDhHmn1T6yOvcDRjdLxD1dvEuayoVrzdxippc3Lz+pq3Bw1VJsdaxqqJWlfoZrna1QJdwAIEOFegYZqnRpoqKaOhmpzaKiOIdwBIECEOwAEiHCPeIMrTGO8wFGvU/P1C3GFao3pqV+h2uYFjV/wUcvPUrMNfhqqJalf1Jzj4SbcASBAhHsFGqp1aqChmjoaqs2hoTqGcAeAABHuABAgwj1hjf7kbxHkp9a81JF/+XnPkBeEeyT9EzP4MKI2snpi+IJswj01cRunQCtYv1CJcAeAABHuABAgwj1hjf7kbxHkpdS81FEERVq/kA7CPcJly8gTwnpi+IJswj01NLzQTqxfqES4A0CACHcACBDhnrBax/o4hto8Riw+1i9UItxHpfvhKOJHkfxID0M9MaldoZrjDwXhnhIaXmgn1i9UihXuZjbPzDaZWa+ZXVPl8ePM7K7o8cfMbHbShQIA4msY7mY2SdIySfMlzZW02MzmVsx2paTd7v4WSTdK+mrShQIA4ouz5X6BpF533+zuhyTdKWlhxTwLJX0/un2PpIuN/UQAyIw1agiY2YclzXP3f4zu/52kd7r70rJ5nozm6YvuPxfN83Kt1+3q6vKenp6mC16xaou+89DmWPM+u+NVSVLn9KkN5z0wdFh9uw/Enr/WsmacdoJOnDJp3OO79w/p5VcPjr5+s8sbef3pJx+nU084tun6Kl+n2jJHHpt95ok6dtL4//dfOzis/r2DR03rnD61qXFupsbysUxyGbVea2T65GNMc6adNOHlTMT+Q4f10p7468eW3fs1OHQk9vxFFmcdnnXGCTp+8vjPYdI11Pq8N/LZizv1/ne8saVlm9lqd+9qNN/kOK9VZVrl/whx5pGZLZG0RJLOOeecGIse77QTj1XnWfFW3u37BnXy8fHn79t9QBe95cyWwvPs007Qg88M6B2zTq05T/cT23TJ26ZryuRjRpf3rnPP1OknNV7eOWecqPuf3qGu2ac3XVu5rXsHdcZJU6qOyUnHTdbaLXs0942n1Hx+/xPb9Be/P11Pb92n/UOH1XnWVG3bN6hTT4g/zo2cdcrx+nXvy0eNZd/uAzrrlOMSWcbQ4SPa8crBca/15o6p+vmGbXrf3LNy8V2cL+05oHd3TtPJxzf+mL5l+lTd++Q2Xfb2szTpmBwU30a79w/JTFXXhWlTj9Mjm3fqD2fU/hwm4Q2nHq+Hnn257ue9nolsoMUVJ9z7JM0quz9TUn+NefrMbLKkUyXtqnwhd18uablU2nJvpeBL3/4GXfr2N7TyVAB43YhzzH2VpE4zm2NmUyQtkrSyYp6Vkq6Ibn9Y0v96nk8ABYDANdxyd/dhM1sq6T5JkyTd4u4bzOx6ST3uvlLSdyXdama9Km2xL2pn0QCA+uIclpG7d0vqrph2XdntQUkfSbY0AECruEIVAAJEuANAgAh3AAgQ4Q4AASLcASBADf/8QNsWbDYg6cUWnz5NUs0/bZAh6mpOXuuS8lsbdTUnxLre5O4djWbKLNwnwsx64vxthbRRV3PyWpeU39qoqzmv57o4LAMAASLcASBARQ335VkXUAN1NSevdUn5rY26mvO6rauQx9wBAPUVdcsdAFBH4cK90Zd1t3nZs8zsATPbaGYbzOxz0fSvmNlLZrY2+rm87DlfjmrdZGaXtbG2F8zsiWj5PdG0M8zsl2b2bPTv6dF0M7NvRHWtN7Pz21TT75WNyVoz22dmn89ivMzsFjPbEX1r2Mi0psfHzK6I5n/WzK6otqwE6vqamT0dLfvHZnZaNH22mR0oG7eby57zJ9H73xvVPqFv7KhRV9PvW9Kf1xp13VVW0wtmtjaanuZ41cqG7NYxdy/Mj0p/cvg5SedKmiJpnaS5KS7/bEnnR7dPlvSMSl8a/hVJX6gy/9yoxuMkzYlqn9Sm2l6QNK1i2g2SroluXyPpq9HtyyXdq9I3aF0o6bGU3rttkt6UxXhJeo+k8yU92er4SDpD0ubo39Oj26e3oa5LJU2Obn+1rK7Z5fNVvM5vJb0rqvleSfPbUFdT71s7Pq/V6qp4/D8kXZfBeNXKhszWsaJtucf5su62cfet7v54dPsVSRslzajzlIWS7nT3g+7+vKRelX6HtJR/cfn3Jf1V2fQfeMmjkk4zs7PbXMvFkp5z93oXrrVtvNz9QY3/drBmx+cySb90913uvlvSLyXNS7oud/+Fuw9Hdx9V6dvPaopqO8XdH/FSQvyg7HdJrK46ar1viX9e69UVbX1/VNId9V6jTeNVKxsyW8eKFu4zJG0pu9+n+uHaNmY2W9J5kh6LJi2Ndq9uGdn1Urr1uqRfmNlqK31XrSSd5e5bpdLKJ2l6BnWNWKSjP3RZj5fU/PhkMW7/oNIW3og5ZrbGzP7PzN4dTZsR1ZJGXc28b2mP17slbXf3Z8umpT5eFdmQ2TpWtHCP9UXcbS/CbKqkH0r6vLvvk3STpDdL+mNJW1XaNZTSrfcidz9f0nxJnzaz99SZN9VxtNLXMy6QdHc0KQ/jVU+tOtIet2slDUu6LZq0VdI57n6epKsk3W5mp6RYV7PvW9rv52IdvQGR+nhVyYaas9aoIbHaihbucb6su63M7FiV3rzb3P1HkuTu2939sLsfkfQdjR1KSK1ed++P/t0h6cdRDdtHDrdE/+5Iu67IfEmPu/v2qMbMxyvS7PikVl/USPtLSX8THTpQdNhjZ3R7tUrHs98a1VV+6KYtdbXwvqU5XpMlfVDSXWX1pjpe1bJBGa5jRQv3OF/W3TbRMb3vStro7v9ZNr38ePUHJI108ldKWmRmx5nZHEmdKjVykq7rJDM7eeS2Sg25J3X0F5dfIeknZXV9LOrYXyhp78iuY5sctUWV9XiVaXZ87pN0qZmdHh2SuDSaligzmyfpS5IWuPv+sukdZjYpun2uSuOzOartFTO7MFpHP1b2uyRZV7PvW5qf10skPe3uo4db0hyvWtmgLNexiXSIs/hRqcv8jEr/C1+b8rL/VKVdpPWS1kY/l0u6VdIT0fSVks4ue861Ua2bNMGOfJ26zlXpTIR1kjaMjIukMyXdL+nZ6N8zoukmaVlU1xOSuto4ZidK2inp1LJpqY+XSv+5bJU0pNLW0ZWtjI9Kx8B7o5+Pt6muXpWOu46sYzdH834oen/XSXpc0vvLXqdLpbB9TtI3FV2gmHBdTb9vSX9eq9UVTf+epE9WzJvmeNXKhszWMa5QBYAAFe2wDAAgBsIdAAJEuANAgAh3AAgQ4Q4AASLcASBAhDsABIhwB4AA/T/9zv9a4nbUwQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(range(n),x_true)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x11aee9630>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzsvXeYJNdd7v85VZ0m593ZvMrBkpAlIWw52zgDBhtM+l3ANvheojGX52IwXNKPCzZgG4yvbTDIlgEbA8bZVrKClbNWWq20eXdmd3Lq6Vjp3D9Onarqnu6e7km7s1vv8+iZ1Ux31amqc956z/v9nu8RUkpixIgRI8a5D+NMNyBGjBgxYmwMYsKPESNGjPMEMeHHiBEjxnmCmPBjxIgR4zxBTPgxYsSIcZ4gJvwYMWLEOE8QE36MGDFinCeICT9GjBgxzhPEhB8jRowY5wkSZ7oBUQwODsq9e/ee6WbEiBEjxqbC448/Pi2lHFruc2cV4e/du5fHHnvsTDcjRowYMTYVhBAnmvlcbOnEiBEjxnmCmPBjxIgR4zxBTPgxYsSIcZ4gJvwYMWLEOE8QE36MGDFinCeICT9GjBgxzhPEhB8jRowY5wliwl8BHj8xx4Gx7JluRowYMWK0hJjwV4A//vp+PnL7wTPdjBgxYsRoCTHhrwCW42E53pluRowYMWK0hJjwVwDXk3hSnulmxIgRI0ZLWBPCF0L8kxBiUgjxbOR3/UKI24UQh/yffWtxrrMBrpS4Xkz4MWLE2FxYK4X/WeBNVb/7AHCnlPIS4E7//88JeF5M+DFixNh8WBPCl1LeC8xW/fptwOf8f38O+NG1ONfZAFfGlk6MGDE2H9bTw98qpRwD8H9uWcdzbSg8j1jhx4gRY9PhjAdthRDvFUI8JoR4bGpq6kw3pymooO2ZbsWZhZSS3/mPfTw1Mn+mmxIjRowmsZ6EPyGE2Abg/5ys9SEp5d9LKW+QUt4wNLTshi1nBWJLB4q2y789NsL9h6fPdFNixIjRJNaT8L8G/Lz/758HvrqO59pQnK1B24MTi3z1qVMbci59/WfjfYgRI0ZtrFVa5heAB4HLhBCjQoj3AH8BvF4IcQh4vf//5wTO1rTMf334JL//lWeX/+AawPPXnTln4X2IESNGbazJnrZSyp+u86fXrcXxzzacrQuvHM/DdjdmBbArtcKPVxzHiLFZcMaDtpsRZ6ul40mw3Y1pV2jpbMjpYsSIsQaICX8FUEHbM92KpdAvIm8DGidjhR8jxqZDTPgrwNmah6/bZG8ACWtLJ/bwNw8WCjbffX7iTDcjxhlETPgrwNkatNVN2ghbJ87S2Xz4rydHec/nHiNXds50U2KcIcSEvwKcrUFb3SZ7A0o3x1k6mw8lx0PKjekfMc5OxITfIryzWNmeCUvH3aAgcYzVI5iVnYViJcbG4LwnfNeT/OxnHuK+Q82tGNWD5WxU+LptG2npxAp/80Cn7J6NYiXGxuC8J/yS7XL/4Rn2nWquJszZ7F3rzBlnA3Il4yydzYezue+eDfj2M2P8zy89faabsa447wlfK9RmrQlPnr2DJrB0NoDw4yydzQcnJvyGuO/wNN95duxMN2Ndcd4Tfqu+pv782ThmNM9bTpylE2Mp4mfWGAXLPSvH9VrivCd8x2vN19QOxtk4aLxAdcdZOjGWwnE3T9DWcT1Ktruh5yxYzqa4N6vBeU/4raqeIDvlLOwYQVrmBlo6Z+OLL0ZtuC2KmzOJj95xkJ/8+4c29JwFy92QVepnEuc94Qeqp1nC15bOGneMp0fmefv/vX9Vqka3bSMsHS/28DcdNpOHf2quyOhsYUPPWbTcs1LIrSXOe8JvVeF766Tw95/O8sTJeWby1oqPsbGWTpyls9mwmTx825OUN3iBWMFykTLMQDsXcd4TvlY9zSpVPVjWumOsxUImzb0bYuno+xYvvNo02EwK/0x4+EX/fJvg9qwY5z3hh1k3rRF+9b9XC7kG6ly/NDYkS+csXoAWozY200pbx5U4ntyQNSUaeb/G0GZ4Ia4U5z3ha4JtVeHD2g6ctSjZ4Hmrf2k0fy78c527g+Ncw2ZS+Lbfxo20dYqWVvhn//1ZKc57wm81CBsl+bXkVe2MrIZA3Q3M0jmbF6DFqI3NlKWjlf1G2TpSSgr+uTbD/VkpzkvCL1our/7Lu3jwyEzLHr63TgpfNiDQbMlu6hgbWh45KONw7g6OKI5O5Xj/vz21YVtIrgdazUg7k9BtLW2Qwrdcr2V7dzPivCT8uYLF8ZkChycXV6fw19LSqZPm+NzpLNf+8W0cncotf4wNLK1wNlcNXQ88eHSG/3ryFKNzxTPdlBVjM2XpWH4fLm+Qwtd2DqztzP1sw3lJ+GHNGRkoiZV4+GuZi+8GK3gre9tEtoQnYSJbbrptG1Hv3N3AeMHZAH1P85t485DN5OHrflWyN6Z/FSKEvxmC2ivFeUn40TKxrebVR/ltLQeOV8ci0W1thlg3cjHU+ebha8W5WNq8hL+ZFH5o6WyMwi9Y4XPdDPdnpTgvCT+6UUir1TKjL4aNyNJpJd9dk7C1IUFb9fN8ydLRcZHNvD1gUDdqEyhYO7B0Nl7hxwuvzjFESd5tcRBUWjpr16Z6BGq34MuHls4GBG03kVrUsFexmMc6ByydzfTM9DjYOIUfWzrnLLRatj3ZcuaCt14Kv45F4rawTiB8aWxcWuZmUvh/9LX9/OLnHlvRdwNLZxMT/qby8P1xeSaCtpvh/qwU5yfhB/nIXsuqZ72CtvUIVFsJzSj8jbR0NpNa1BiZK3J8Jr+i7+qgbS728DcEgaWzQWmZhThL59xF1BdvVfVU5OGvA+FXZ+m04uGfCUtnI5e+rxYly12xB68JaDNbOpsqD19bOhuk8KNB2zgP/xyDHUnFbDW9cN2CtnU8fKeVLJ0NTJXUl74ZyEOjaLvky86KgnJ61rSZg7abqZaOnlHFaZlri/OS8KPqVBNssxy5bpZOnZmGbl8zq2c3srTCZtzTtmi72O7Kyu7qgnSbOS2z1d3dziRsT1s6Gx+0PZc3QTkvCT9aMC0MijZHAhsdtA0WhjXl4aufG1JaYRP5wRraHliJSj8XLJ2Nfmb/eN8xHj46s6LvBnn4TSj8Z08tMJktreg8GsVoHn6s8M8thCQa8fCbfMZR3l1bD99vWx2F31SWzkaWVtiECj8g/BWodPscsHQ2Okvn4989xBcfHWn5e1LKljz8997yGJ+850jL54kiDtqew4iSaKh6mvTw1ykPv5760sr+rLN01snaWk/o1LuVkLbOw4+mZZ6aLzK5uDpluZHYaIVftj3GFlqvPRQVEc3Yb3nLrUirXAkKkRdLHLQ9xxAN1IaZC819d70sHVlHMQcvp2YsnRb8/lZwar64JNAZbeZmUPlSymBHoxURfg1L531feJI//vpzNT9fsl3+++cf48QK00DXA0G8aoMIrey4jC+0/kKMZqQ1o/DdiHBbKeI8/HMYlR7+yhX+WnaMcIvDynbottrNWDqBh792Cn8iW+KVH76Lu1+YqjzXOt2H9YLlesH9WZWlE/nuVK5Mtli7dPXoXIFb90/w+Im51hu7TnC92qKiFTx4ZIZ/efjEsp9z/Ps9tlBqOSvKjozFZhS+43mrFl/RF3ns4Z9jiAZCm/U1Hzs+y7HpfIU6WtvyyH7bVqHwwyqga0f4cwUL15OMzBUqzxW59s1QMbNkhW1cjaUT/W6+7NRdH6Efwdm0X4DuQ6t5Qf/Lwyf46O2Hlv2cJuqy4zFfaG4/B42VKPzV2opFO66lc84i6mU2uwvQ+7/0FB//7qF139O2bpZOE+cKPfy1a5dexFU9aNfrPqwXogN6JeUR9D3NW05ALrmyU/dlF87Mzp6X4Vp4+HMFi1x5eQKPKvOxFm2dqGBpJkvH8eSqbcWC5SKE+vcmWkvYMtad8IUQx4UQzwghnhJCrKyQyRpDD0K7IkuncYeZXrQo2966BSvrTbdbUe1yHYK22rteqLIuote+GTz8KOGvxtKRUgX4HNejZHt1X67e2ajw14Lw87Z/3Y37WDR/fjzbWuA2euzl8vA9TyLl6mfbBculM5UAVnZ/siWbP/zqsxUrds9GbJTCf42U8lop5Q0bdL6GqFT42juv/5BLtusv2vHWMQ+/sm0aQT38VkorrJLwpZTceWAC15PBsaoVfrSZm0LhR4JyK8mlj9YnypUc8lbj/U8DhX8WycW12MJvvmABy99DaxUKP9rXlyuP7NaZGbeKouXQmVGEv5L7c/+haT734AmeOjm/qnasN85LS8eJFCRrZrm5VrfqBRH+fl02QKmj8JfzyaWUYRxglapy/+ks7/ncYzx4ZCYgrIWiVdmuCg9/ExB+VOGv0MNPmiL4via8eoTurYO9thpEc9tX87zm/Bf/ciuOo5ZOq5k60b6+XHnkULy1dIolyFsuHemVE/6peTWLOdvXaWwE4UvgNiHE40KI927A+ZZFdIl5M77mnK9qbE+2tKet5Xj8ws2P8OyphWXbFJZWqOy5YbXMxueKNn+11TJ1p81FgpJLFH7Uwz9LSK0RosG/Zsoj/N13D7FvNFRrtuvR154CKgm/Hnm2skJ6IxBt5kqtSD3ThSYI316Nh9+8wnfqjJtWUbRcOtMrt3ROz6trzDewdKSU/MO9R5lY5arg1WAjCP9lUsrrgDcDvyqEeGX0j0KI9wohHhNCPDY1NVX7CGuMoD6N11y1TE12jutVpSM2Ps9s3uLuF6Z44uTyqXn1s3S0pdP4ZNGXz2oVvp6O265X18PfbFk6UUunmaDjR+84xDefGQv+33Yl/R0+4ZecIPBb77kEAfSzZPYTfUYrVfha+MDyStZyIx7+Ci2dznRieYWv19Gs4jZLKSlYDl2rsHROBwq/fnuncmX+7FsHuHX/+MoaugZYd8KXUp72f04C/wXcWPX3v5dS3iClvGFoaGi9mwNEOkmkHn6jQaB9S8etXOCxnBJoZVof1NJxqwm/OeKItmW1vnGU8AMPv0HQdjN4+Jo4ujOJZclKz/yiL07L8ehtTwLVlk69tMyzS+GvRVbVXD7sA8u9NLUyH+hItbzaVidVdKYTy6Zl6hfZahIoyo5aM6AJfyWPTFs6jWIb+p5YG1TjvxbWlfCFEB1CiC79b+ANwLPrec5mYAeDUTbVYQKF73kt5eG3sqirnocflEduQeGvmvDdsGMGHn7BrshP9jabh+8r/KGudEMVBpWb3GtYrhcq/ApLp47CXwO/fC3hrAHhz0cUfrMe/p6Bdiay5ZbOo0sjd2YSyy68WotUU903OlaRpXO6CcLX/epMxnXWW+FvBe4TQjwNPAJ8U0r5nXU+57Jwa6y0bTQw5wLCr1T4yxF+Swrf022rE7RtwcNfbYeqVPjqWJafhhi2a2kbz2ZopTjYmSZXaqxOw4GpUzFVtlJo6djBS6Petevfn0k1F0V05hi1454emefkTKHWV5Zgtg7hf+nRkYDwNELC7yBXdlhc5p5HocdiMwp/LbJ0dFvbUybQ+sKrku0yk1f3ptHs0WpSvK0n1pXwpZRHpZTf5//3Iinln63n+ZpFdPWqJtJG5K2Vje1WBm2Xt3T887Vg6dTd4rBJS8cQa2fpWK6sONZ8JFNn0yl8O6rwG6tT260kEdfP9dZB27zlNm/ptBjfODlT4AP/uW/NSSH6jKKz2d/60lN87I6DTR1jrhC1dNT1l2yX//Wf+/in+45VfFbnz+8ZaAdCy6MZ6D7XlUlQsr2GBBzUwlpFqqnu7xmf8Fs9VvTaGip858zHdc7rtMxmFX69oO16WDr1tzhcxtLxP5dOmKsm/HJE4doRhRrN1Kn0hM8OFdsIRb+0glL4yxF+5dRbK7OOdIKUabBYcgLCqxu0bXJmVo37j0zzxUdHWs5sWQ7R5xXt6wXLbXrl8byvYg0RLl7Tx9o3WpmJplXzRUOdAJyaa57w9T3TFkujrLO1sHR0gDmTMFd0rNMVhF9/RqLPcybXZpyXhB/tJLrDSlnfx9fZCUvz8Bufx2shU6Nelk6zC6+0KskkDWxXrqoeSGDpOJUrSSsIv2Kms+JTbRiKtksqYdDTliRvuQ1jNvr69YtMK7OUadCZSZAr26HCr2fprDAPf71KGEdnGtFrtxyv6X1j5wo2nekEXZlk8MLTVtGzpxcq2qwJ/8KhDqA1ha/bqhdCNSqvsBYVQHVb23yF3+qxNOEPdKQaWzrOmQ/kn5eE70RKK1Qo1ToPWhOdXVWVb7nMgMCmaaHwWT0Pf7maLPpcaV+lrMZmCS0dr0JdRVMzZYWlc2Y68M9+5iFuvv/Y8h9EWQ9tSTPIxGiULx0ofO3D+/+fTBh0phMsliJB2+UUfov3ZqXfa/a46tiVwejmCd+itz1JZzpB1vfk9XgoWC6HJ3Phcf0+tKO3jVTCYLQFha9fkjovvtygfSudSdU6Xyah6LDVW39qvoQh1GymkaVjVc0czwTOT8KPpGVGB1Y9VaW9a8etrMq3nNcXljpo3sNfUlqh2aCtfxnppOGfc+WEESX86Hmjq23PhuJpz57KcnAit/wHUZkYbUkzIJFGSizY5N4N7wNAyhQMdaWZzJYDG8SrMzNcKRGtxWrYRseFyn5ru17FKuRGmCtY9LWn6MokIpZO2M+eHgkXqmkPP5M02dHb1pql44UePqiXSb2XUpBltwYz2rYWPXzPkzzqV9Hd2p2hpz0ZlNyoBTuSDHGmcH4SfmQwNkNcOljVctC2TuZNLejDLlX4zXWSwNLxFf5qVETgNTpVQdsKDz/8/ErI6SO3H+TD33l+xW1UbWi+LG7RdmlLmYFN0MjHr07L1AM1aRoMd2eYyJYqlFyt63eqZgfNIsggW2MVWK/on+V4Te8WNVew6evwCV+/8CKX93RkZbLOOU8nDHb2tTHaUtC2UuH/1W0v8Ja/+V7Nz66Jh6+Dtknf0mnyWF9+8hQ/8akH+frTp9ne20ZnOtGUwj+TBfXOS8KPBmrrKR8NKWW48MprLWgbZt40b+kszcNvTvEFQds1VPh64VUqYZAwRIWlU1FEbgUd+KGjM9x3eLrl7y0U7TCu4XlNv2yKtksmovAbBSqrp976fKmEwZbutE/4IUnWer5eICpaew7rpfCj/SG6utyTzZUgBpWt1udbOrka6xCigVvL9UiZBkIIX+E3l/oZbat+Od/1/GRdS6iZWljNnk8TfrPH+tYzYwx3Z/jpG3fxcy/dQ3vKbC4P/wwmOZxThP9/7z7M5x9afjeeKGFUKPwaxFWw3DBNr8WVtq0E7upZOs1ugBJ6+Ib/+TUI2voefto06G1PVqy2recJNwvH9VZUtfKNH72XWx48EbSh2Qwh5eEboaXTSOE7lQq/XKXw85ZbUQ+l1vNdKXFHV4GvJWopfD0OmvXwZ/PK0unMJIM8fN3Mwc4Uz49ng9hO2faCvrizr43pnNX0eZwqhZ+3XCzXqzne1kLhl6sVfhOHWizZ3Hdomrdes40/f/s1vO3aHRUvwkbniRX+GuGrT57mb+44tOyUrELh11mQoqEzdHrakthVinJ5S6d5ladPXW+l7XIvDd0W3WlXpfCrVtom/eyWhUIdhb8SwvckhRY3npZSMp4tMemTbSsbXxStKkunCQ+/egFWyjQY7skAcDyyV22t5xsuuluZwl/rwJ5T4wWtCagZInZcj8WSEwRtF6s8/P6OFLYrg75Tdtxgtrmjrw2g6cBtoPB9wteoVRt/LWr86zaHQdvlj3XXC1NYrsebrhoOfteRViuD64336v50JnBOEX7RdpnOlXlypHGxsmgqZtRjrdVptG892JlSQVspgzK5y1s6ledrhHDFYGVnCJViawp/NRUzy1VB24Qh6GlLLll4ZRqiqbbVgu3KlhV+4IH63r2UzQ/0kqOCtm3+C7ERyS3x8H3yTSUMtnZnKn4X/VwUYbXMFhX+OqVl1lohrq+zaLvLpvHq2V1/R8qvR2RXHKvdz5kP+o6jLB2AnX2tLb6KrrSNopb1FNT4XwMPPwjaNjjWfMHi43ce4u/vPcJgZ5rrdvcFf9Pllevl4ocz51jhrwn0IL51/0TDz0XfwFHV0Ijwh7rSwUIt3ZGX49RWMjXCFM46Hn6TpRXSPqGtpaWTNA1621MVHr7rhS++lZCT63nkreWJpla7XE8GPmiz5y5aysNPmsvHOAIPv8r6SJoh4UO4FL9WLn4razCicIO4zxor/BovKH1dnlyehPQ46GlL+iUP1OxPt1PfCx2sLTte0Bd39CqF32ymjh6fHUsIv4HCXwMPvy25fB7+bc9N8Ne3H+TgeI6f/YHdgegB6Eyr7+fqpPxGreQzhXOU8McbL8eODKay3Vjha0tnsDMN+Bth6KnfGgZt9anrefjLZunooG1i9UHbckSJ2K4klTDoa08xkS1H9t4lePGtzMNXL8/limNFEbVaWlXCJdujLWmS0jOgBucNF7t5FZ9NmoLhCOH3tiUrPhdF8NxarKWzXlU29XFNQ4QZRJG2LZeaqWdjnelERaaTGxC+Vvhu8FP3xa3dGRKGYLTJwG11lo5GLcJvdk/qRliSpdNgXOs1AQ/87mt5/+svrfhbqPBrE74Vp2WuLUq2KnB1YqbAkal83c/VWhFY/XsN7VUO+IWztOKt9/la52mGEL06nw3q4Td5Lr3wajWWTkWWjr/T07W7ephaLHPCL7TlSUlqhUvRIcxUaMXHjyr8VoOiOktHPzurgaKttnSsiMJvS5l0+4TX49fWqaWOV7qAqpVN61s6rhemSVYHbaHx4iYISawjnahYyxASvq/wnVDh65eraQi2dmeaLhfheB4JQwTHvHy4C6ht6YRrataO8BsNHd1vdD+KomOZNR5WIFhiS2fVcD0VMLp0q6rdMbVYvyRrdHobVQ21Bpn+u1Y1ZdsjaTRnZdSzaRp9dkkeftOWTqXCXwtLp+yohWlJ0+BlFw8CqtaLbmdwrhUqfGhtf1ndLseTkWyW1oK26WYUvlNJiPqn/q4O3AYKvwapu67Hjxj3Yzit1cRZC8Va+7hhHCJU+OE5llP4uYjC1wuiFiMKvyNdaelYTpilA8ofX25Dcg3HlST8RW6//9YrePfLLgBqb3e4NrV0NOEvP3PX/SZVg/A7m1T4cWmFNYDuTL1tSnU12iknOkCXU/hh6dRwymqaAiFaKY+88jx8rYSXL62gftbKwz85U+BT9xxp2i+3IkRnuZKEaXDBYAfbejLc7+fOR4PXbo3rOzlT4Pf+65kGGQs+4TcocVCvXU7EO25uUZtcovCb8fCrSxzr72ofX2+IUuvl2p0/zt+mPsHL7PuXv7AImrXwWoU+bso0gn4ZnQUul4uvn1OHX0sHVGqivkdtyWpLxwtmmwAJQzStbC3XI2moHP5ffMWF7Oxv89tYg/D9a1nN+3GpwpdkS3bNrUmjazKqoYu91SP886Ee/oZBrxbUg7DRXphRkij5RbWqf6+hO3BQ18PxMIXAFKIJS0f9bEYBhyttq7J0/M6xXEZKkJZZw9L598dH+ItvP890zqr53WpUWzopUyCE4KaLBnnwyAyeJ31Lp77Cv/vgJP/68Mm603h9nY2qC9ZrV7TKaTOEHxTHSpqYhsA0RFMefnV6ZjJRh/BrtCFpZwHoceeX/K0R1jtLJ500lgRtoRmF728SkjZrWjqBwg9mh24FKSZNo2llqxW+hibiWmN6LWoPKZtWkDDC7Lt/eegk7/jkA0uyf2zXwxBUBGs1wvuyXJZOrPBXjZJ/M3v8aXaj6WPFJsmRqWc9hZ80RUWwzzAEhiEqMwOspTGDejXua6HWZ6VUXrVW0o06iv5+pkaWznHfd5/ONbfzUFgtU5VW0Mr25ZcMMFeweW4s62fp1L9vOvZR7zno9hVaUPhRbz2MbSw/eLQYaPNnP0lTNLyXSz38UB0DQeC2x59N1iKyhKNq/HR52WXbF0W07MdaIqrwq0tGwPK5+DWDthHCb6vy8KstnYQpmrb+HM8jEbFMtIipmYev97NYBYdajurjQvgzd0+yWLIpO94Sp8CKxCaqoV96yyn8M7l/xLlD+H6H7WlR4UPozdZK7VIrBs3g7R9V+MHbf/Yo/PkuOPVExXdbWXhVa+eeatXeqKM0ytI54S8SmmlW4UcsnSjh37CnH1BL6D0vnNbWum8h4dexdAKFvwJLx5MtBeu0etWklDKNhtlBQfG0oDxypW+7tadS4deaoicddc+7ZWuE30qgv7Xj6jiEGa4gbkHh58sOhqCiAN1iyQmevbYzdPA3KqTwlEXT7O5ftiuDOBmE3nqjPPzVpmXqvmwIJeT0cavPGU3aqMayQdtY4a8dAkvHV12NFL7jeYjIjCwdZJssfRA6vUxPMcuOG9gCwXObOwHShckDFd91WwjaBittI5/Vg17vxJMvO/zx1/cHpWmjCIK2NTz8EytU+GW/Hr7u4KGt5eLKyHqEGtenF+bUI9YwaLsCSyfi4TdDjJrM9OwnlTAae/hVC2TC8siqD1y/u48LBjuCzT1qzTK0wu9ZIeGvdWkFfb9TidoKf7ksnVzZoSOVQAgR9Ifoc6hW+GXbJ9HDd8Jf7KbHyDWv8N0w9RnC59YwD99b+R4Q0UViyqoNj1t9zuhnq5H2a07VDdrGxdPWDkHQ1lddy22akIkElEIPv9ZxPf9BGsH/G0JgRIO25UX1M1e54Ev37+by8Jcq1oDwfRJ/4sQcN99/nEeOzi79vn+K6mqZ8wUrWDDVLOGXI0rE9v1NCD1sveI46uG/+7OP8u1nxoJjBAq/xnPQVhW0ZulEFX4rXree2eiFNUmzsdqsVy1TD/Qrt3dz12+/msFObenUV/h9rIzw1zqwF83SqVXjqWi7fP6hE3z1qVM1v18ou4GCDVdYh9VKdQpltLR2OmHCzGGwFtnujjft4dueDGbUUEn4UsqK40RfjCudFEVtGsNQY1Gfo3rmE50NVEMIQUeDipmxwl9D6C3sepvw8B1XBkoYIqmMNRW+WjGoSa9sexGFrwnfH9S5yYrv1sutr4VaU3mtnDVRaWVfK8feraPwj0c2qG4+aBtuxRa1dPQgtP2qofo+lB2P7z4/yb2HwuqXuYiH/8zoAtf/6e1BqmyUpBvVD1/arqUe/nKEf2Qqx2984Ul62pJeDb5jAAAgAElEQVRcvbMHWF7hV6dj2q6aEVYH6pINFp6lfMLv9RX+l58YDaquNkKz19UqdBvTEYWvy2CDEkife+A4X36iNuHnLIf2tH5hhmnJS1baBgrfX3hlqZnOIHNNv8ScKtsksHQcjy8+OsLLP3RXoOZbqW1VD2W3UuF7keuqVvjRGW8tdKYTdft0XEtnDRHNlzfEcgrfq6nwawlx3XG1wrdcjx3eabaLmdA3rKPwWymtUCtLR/vcWuFki4pEa6nTJTte+Z1K+/dCtGDpRLJUoh08SGl01L4AhlCEr19Ek5EKkrr8cNnxODy1yEze4tCkuk/RQbqiPHy3eYX/f755gJLj8sX3voRtPW3BdTRamFa9iXnZDYN6UQRKt8axAoUvFpnMlvitLz3NN/aNLflcNdbPw19K+HY0D99ymS/YdWej+bITWHq1Ff7StMxUwgBLCY4+b67pTBq7OksnESr8I5M5xrOlJc8IVr4Jiu204OE3CNoCDUsknw2WTmL5j2wOlCI77GSSjRd5uJ4M6l5AEwo/YWAGCt/lfy7+JTNeG1/x/s7/kCb8KoXfQh5+rSyd6gqY2pqpSfhVQdubHzjO1/eNceNeFWi9eKizZQ9fb3GY8r1r01BWliorrQZHwhBk/XZNLEYIPxK01YNGK/zo/WjFw69VB385YpzJW1y9o4crtnUHv0uZRsWio2pErSMpJbYja/q2QfZUjTakXV/hk+dwXt2XpqpSthDobwWhwjcjCr8yaDtfsCpeAlHkfQ8fCMRPPYWvf59OmEH2Wr87o4hucRy6hmucIYTtesE5AAxDkDINSrYXBEQt31qpleTQKiw3aukoha9fKEs8/AZBW1CB23pB22BBX1xLZ/UIU+/UasplPfxkDYVfK0vHcUknTJIRD3/Qm+IqeQjp+Z0hsHSqPfzmfeZaVf+qN2bQSrpcgwyql7gfm87zyLFZbnnwONt6Muzoa2uK8B3Xw5NqRhCUR4508ISpNklXlg6+wlcdfHwhPH4QtLXD7el0DfmowmnJw49YOs0qfL2XbRTJRGOFX+kRy7q+rSalWuSc8gnfEJJCdgaoH8COYr2zdCqCtpF2z+QsHE/WJaNcxMPXzpbjyUi1zDBXXj+ndDK0dHq9Wb7PegI+cgXMn2zYVscNU5E10kmDku0GZKqDzLXKPreKiqCtoRW+uobaQdulOfgajXa9KkdmzgEOfB0++0PhFH+dcc4QfinSyZZT+MrDjyr8+hUmlYcfZum4nku3t0AnBXpLvt9Zqk34ejw1p/D9ttVQLDp/vKHCj5Sp/cTPXMc3fv3lXL2jh2zJYc9AO4Od6abSMjURdvpqrmi5FWoraagcdl0eOarwZ/Ll4FprKfzJbHnJNTaqS1+vbbbbfFpmyd/aMIq0aTQsahYdkI4n/TztpYNc94la/UYrfAA7q2Z+zaQlhrV01kfhp8wwaBttj34Z1+urytJR91H4MzvHDfc8TpoGSVNQdtxg7KVMI1D4ve4Ml7iHQXqwUDtO8I19p/nH+475tXQqqUmP6ajChypLZ4WEHxU1hhB4Mpql4y35bCNLpyNt1p21RjfWCTKKDt8Jx78HdvM7gq0G5wzhlyOpd8srfC/Y7ADC7IuaCt/fuUcP+F5ymKhj7yz6aZja0ilnA88yerzqxVS1BpWsMRvQxBN6+Hpv3VqEr36ahuCt12zjRdt7+JVXXwTA3oGOgPCXS13TJNARWVms0xFBKXzH9SIevhG0S0pl20gpI0FbL1T4i5rww/avvHhaaLsEmD0Kn/uR8AWMXzQtUa3wGy+8iqp/x5NLZjnBcfzfifLS1bQpr4AnfctncXrJcethh3WEf0j+NdJpzn5rFjoBIJ00gn/re9CVTjC2oEoXawI/NLHIPQengv4oSgv8wvj/gdljAEHSgn55mIYgnTApO14wk0knjYDIup1ZdnqnVWPKtTOXvvLkKT77wDEVN0pUE75v6ZQq41gNtyj9xvvh9j9c9t5UZOn4C6/0faiVpbNSSyf6/ANRoWc7paVlHNYD5wzhRy2dZjz8dA1Lp9aUUFs6WnEMiLCz7qomfFAq/4lboJyrSfhfeGSEV374riXEG62DXk3+YZZO/aCt/mw0keQNLxrm7dft4C1Xb2OwM4XlekHgtx5Cwo/cn0gHT5oGth+sCzz8yLqAiWzJL7rmBzwdN4ivTNawdFa28CpaSydyL0YegWP3wPTB4Fe6aFoUqeWCtpH76wRxjFqWjuDF4hA/fscrYOpgxd8ybp7TDADg5RXhN1oMqPGqwu283nycrkJtFbxS1FT4ugxxJsG4XwJD35c//eYBfv6fHuH1H7mHycUSP2Z/k2vnb4MTqjZQ0jQq0mMV4RuUHTe0dCIefpczww45rhpTqk34ZcdjNmepgn1VGVGZhFlh6YQv/0haZvX4PXwHPP5ZcBv3serKntHYRC1LZ/ksncYrbX9AHIC7/1z9Mib8laHkL4hKmkao8HOT8A+vhbnjZEs2b/joPTw3OodwrUqFv0xphXTCCDITBn3Ct0mwq/iC/6FFwO+gz30Vvvbr8MyXIlk6YaccnSswtlAKyFtDei5bUfn1uhmqg0j22EpVNWPpGJGBYhqCj7zzWl556VBQz38631g5anUWrUWerCB8ESh8nZYZfYlMZEsVL4Cy7QVEN1kjaLtShR9UEY0+M/3iLYaKu2R7FfEafT3ReziXtyoKZdlVSsx2ay+2SRgGP2AcQCBhodKXTnsFRrwtAMi88vCjaZD1cIX1jPq+3VoNnuXgepKkIel2ZyuKwqVMg7akycRimS4KJByl9AtltZ3h0ek89z93gp83vqUO5BN4oPCrCd/2QksnkpbZZc+yR/oKv1T72ixHbYqTL7sVWTqgZrkdxVOUSqp95RoKv6IvSKnGf2keTj3e8N5YkecbZunU8fBduYyls3we/tvN75G678NQzsHCiPpjTPitoWSHNk1aK/xTT6iHffpJxuZLHJzIIe76E/4j9UcVJNColo5lO/SJbEB6g6gHs898EbvKh5R6KGehd5f6wnNfVT/HnwkDsTJUH/p3c/lKP/118mHuTf8mgyxU5GK/xniKX37+57hGHAmsk5p5+HrgidoBpYDwG5SNjh5b10sBKgZfwhR+0NbP0jFFRXsmsuWKDcKjlo5W+FGSaEXhR4ua1ayWqQeNTyi6ZHZ10DZVFbT94Fee4RdufiRynso4Sj1VlzAFLzKO+xeaq/hbxstzUirCF0Wf8P0Bf2QqVztjpzjPBc5R1cY1JnzHk/yUeRf/4+l3kPYUaepFdZmkytz5bOpD/Hb5E8HfLt2q6tD3PXMzfcK/Pv+lmjCEn60V9js17sKYTTqSlmngMij851PH0tHPZCJbqqilA9CdsPj/T/8iP1T+ZsVno6re9aSKDxTnVF/QpakP397w3kR9edMQyAoPv1rhu3VX2oISSrYrazoMuv9uE6o/cPx7YRtjwm8NxUhwLlD4WX9aXAxzgNuyx7jGOEaPEaYQNlL4r3Ae5HcPvINMQakTbek8krqRlCzD9AtqEAxcor5w2q+nM/5MReDddmz42m8wtLgfUOmCAPiZPrvlGGnhcJE4XVEB8AcNdbyLxallFL76adQh/AF/Vehyi68CSycVEv4SSycI2uoXjOTHjO+xU0wxkS1VeJhlJ8zSyVtqSq4JtTsTmf46Zbj1gyptb5m21c3S0QrfJ3x93rZUZTdPmeHCq8lsiVv3T1TMuCq9Vl1ewr+vn3873PNhwCd8cdz/UmXxvDavwAzd5GSGRGk2aH/ZcXnr336Pf3t0ZOkFnnwoiA9l7LUlANfzuMY4SsorMejPJG3X41JzjG3mAgMscL1xiD2eapflSnrakrwz8wivHP00t7nX44mEuk67xD97v8NFCw8F9z8hBO2m5wdto4Sfh1RXZWPqWDr6+RYsl7RwK+Jh28UsGVnmYvdIxWejqt6TEv71J1U/ChIoBNbztza8N0uydCIe/tKg7TIKPyiDspTwy45HJmmwVfh7bh+MtCsm/NZQst0g2ybjKw2y/hSyOB963JYihd3uieC76QY7N+1xR0hIi64TdwIwIBbwMDiYuEx9YP6kUiz9F4CI3M6J/XheSCLy5EPwxOe4bPYewFf4j34G/uoSsEv0+TOH3cZEmIvteLzKfBqAnWK6YjehamilY9R5olrhzzSydHKTwUCqa+kYRrDwSVs6P2veyUdTn+TDmc8yni0FGTq6rdFBM5ktBS/f7rYkBT0wjt4DD/6dSlOrg+ggr7kitcrSqa6hE70efawvPTYSqHgdO4laOnqWkEoYqmbSkTthTD2TpJ3jQsN/QVkRhe+USeCQkxnmZBfJshrglutRtFze4t5NeeakIrRPvAT++cdh/Fk4cR+uPyTXmvAdT3Kh8EWLXMDzr/mj8sO8f/EvuclQQmQo8jLY6Y7y5/wtzxiX8z77V3ES7YrAc+NcIY/wxtOfDPZCSD99M59ffA+WbYebhGhLZ+DCysbUUfjRfv3D0/8IN785+H9td17A6YrPLsnDz43D1POBcFjc8XJSk/s4cuxo3XsTDdoKQdXCq1pB2ypRdfopmHweqNrm8OBt8MTnK77bkUqwTfilUQ7dFh4jJvzWULa9YAm28hLdCOHPBSSasBUp7LTCDlCv6qOUkkGppl/tJ9S0cJAFFs0e5swh9aHFMUU0mV5oV7tCsfcVYBfozIe+rnj+GwB0Wkp5LCwuwN0fgsIMLI7R7xP+TjEV+NPJ+SPsFCrgt0uEi7pqBm0j2RK10N+RUqtt61k6I4/CX12KOfUsULmBdGUefpiWKYTgMu8If5j4HFnZzk3yScTs0UrCtytLzE4ulrFdyfsT/8FvereQtxxFtMfUi5DZBgMzkrIYzVcPAuCBwl8AKXFPPATIJYSvSiuoQf2FR0KlXV3/Xp8rsHQOfsdviCJ3Y3J/eNBo4N7/d442ZukibfmE73jYc6f4SOpTvOmFP1Av/KkDcPIh+NTL4KFPsd+4HEuatDkrIADXgdv/N8wcWfonTwZkOSgWcKXEcly2Mc1V5af4GfO7APSTBUel17648AAmHr9c/GWKZPCSHera/evbXj7Czmn13BKjD9PvzdJengrsjCBoO3Bx0A7ZsXVZhQ+w1R6BqReC/PStUo2Di8RpQNbO0vGkstbmTwYK/8SW1wFQPBV5VjXOq/u4Lq2g1yNUZ+nUtPe+/j647YMAFXsF8Min4d4PV3y3P2nRLfyZSzYSmK8T11hrnDOEH7V0QoWv8+RDhZ/yqxhuK4WDop6Hb7keW/zpV3rkATKUGRRZcmYf2UQfHkINLulBphs6t6ov3vhLAAws+kFdJIkXlPfYbU0BMHzo3yDvk3huQg00YLeYDDpx9+jdAOTbtgXEr9tVDR20refhm4agpy3JbL16LmNPAZLMhFKvUQ8/aVamZdq+pWIKwZvt27FI8m7zz3BI8NLZrwSWjipDrCydfn9P4IlsCcf1eL3xOG8tfp12WVAzgGP3qhNUEb7jerzxo/dy6/7x0NKJePjq2v1/aOVYmofj97H133+E68ShpQuvfIU/uVji1HyRPQPtQCTP35FkKCPwgrTMlGnAC9/2z+OreV/pqy/nlS11+x/CwigAOdnGnOwKyLvseEEJ7d25p+DOP4GLXgu/uQ/e8Gew+yV8JfFG5ulamcI/8l24/2/gvo8u+VPSmg/62KBYUM/PzpNG9YeXms9hS/8+LY5hOx5X5h9hLHMRp6XKNpJVhO9hcO2JmwEwZg8B0FseC9R3Rjjg2dCvFP6k7MXrGq7v4UcIv8PNglMMrLIhT/X/TlFimNmKvqDhurb6Tn4qSB+dSu9Wba2KseTKDm/7u/s4MJatsGlMQ+BVl1bIT8Oz/wlS1s7Yyk9BVpXNqFD4hRk105ASd/w5XimeZHdCPVfPVDNu2gcg2R4r/FZRiuRbKw/frfTwdU67ozrrlsLh4LsVhJ+fBjcs7btVzFFOdCLcEjcZ+xkQWXKJXjASLBo9qhogsCjb1JLxzq1w6ZvASDKQU6l6V4ljmIujkMjQbU8CkqtP3gLdO1QDFscZ0JaOmAw6W+/Y9zjibWN24AZ2GaHCr7VoSFs61fVeokgnjLpL5zXRpufVwK1n6aR0lo5v6VzkHmO/3MtEZi8H+l7DG6zbKS2qKetAZwrLt3R29ytSnVos43qSITFPCpvXGU9QWJiE8Wcq2qGRz85w3fRXOH78aM16+BB5UUctHT/d7WLjVN2grU7lfVnyEL9kfiMkEcfmnvRv8W7zOzh+lk63UYTj96kDaL9+bB/Tsoei2aWI8NTjcP/HYN+XVNtFGzN0M2BPYKBmCsb4kzjS4GTmckWGr/kgtPfDTb8Gv/ANbjVeyZzspN1dAQHs+6L6+dxXwS5W/GmgGFqYmvAzZSU+bKFexnd51wIgs6dJOjkuKOzj9ODLgu/JVKd62fkK/ZmOl7Itu48ecgh/HPQ74Ys5gx8na+unnOjmmBzGS3XXJbeokAmuP6/aOOCF/f8i43SQ8VSxxiVK6qOPQKKNGdHvHzyn4mUPfgLsIidnCjw9usCTJ+crSyvo8sjR0gpP3AL/8W44cmft8sjFuWBGUVETvzADrgWFWeS9f8nHk3/HzoQaG/ltL1Hf7d0NmZ6Y8FtFyd+zFMJVeUs9fEnGU4N1MH8Y8Bej+A9bOmX4+PVBjmzZ9hgWc4xueQ0y1cHrjCcZZIFcoh/TEMwZA0HO9we/fYKDV/4G/NinIZGGocsZzCmF/3rzCaQw4UVvp8eeYoh5uq1JePH/p867OE6/8GMLYjLwp9sXj/GsvACrayfDzGKiOnnDLJ06lg6EAdea8G2AjE/4HSmTnzbv5P2J/6i0dAx1DCnBxGOvc5TnvD10pBI8f+G76KbABQf/EVCEr7N0hrrSZJIGE9kStuMEavOt5sN4R+9Tz2L3S2HuOJ7jcPcLk5Se/QZdn7yWP0/+I5ef/nLwolOWTqXPrh5YROH7s6fdYrJGHr7a4lCnhP5U+d/5QOILWH7K3zZnhK1ijkvEaBC0vab4qCLovr2hXz91gINiD5bR7hOhP2j9NEDL7OQ293p65TyvNx5XxDLxNAflLj617U/hp74AO29Y8hzn6aTdba2sMqUsPP9N2Hq1ug8vfKviz5rwPQyGUJZOh63I54Ghd5KV7fyr+1r1mflTXOvuw8Rlceerg2OIdKd62fkv1v2Z6wB4g/kYwl9gNWiP4xXn+UHjcTKeT/ipDk4PvZy73WvxUl1NWTrtjv8Zfw1Dnz3FmFTkfZE4XTMPX5QjgfPRR6FrK3OOeplRzqnZ1a2/B4duC2ahM365EV0uISiPHE3LnDuujnHnn+C4TqXCt0tqcVlhBlw7spG5CwXfq18cg/mTdIki3+eqtTvZ7a9Qf4sJf2UoRvKt0wmDlJ0NlysX59XqWiwSuIx4Q6TcHDtQnUk/wP6F5xRZPP1v4HmUyyUGyFLq2Im799W8xnySAZEln+jDEIJZsz/oDIuyjROZy+Ci16hzbr2SgYKaVu4SkzhdO2D4KlKyzLWGbyftvBGMBCyO008WiwRDYgG7mOexYzOki1NMyl7cnt0khBekc9X28NXPepZOcJ12XinQ574WdMiy4yJ9Zd2RVW3rSCd4i/EwP27eU2XpqLRMV0qGrBEylNkv99KRTjBwyQ18zX0p15/+AruSC3SkErSXJ0nYi2SSJlu6MkwuljEKU5hCUjK7eJWxj47HPwHJDrj6x8G1+J2bv80v3PwouTs+hNs2xIJsp600ieV6fCDxr1wnn6vKv/bvhyaS0gLklDLcLSaD2E7FfUCVf0hjcXnpaUwhcX0b4GJHvfSGfCVsOR4vy34TuncqC0YTfnGOBdFNyWhTv9ODduwp9ZzMdm73bmCMId6V+A5ly6Vteh9PexcyLvvg8rcseUaOJ5mTXcrSaAUHvqZS/N7616qdT3+x4s+D5RNYJFjo2KsUvhsS/gtDb+Sa8md4wrsUAC97mpu8Jykb7aQueGlwDCPT6Vs6qm3PpdWM4B3mfeF5nAkuOHwLn0n9NW1FvzpoqoOHXvwhPun+CG6qqwlLR9KmCb+gxmiPPck+70Kysr2C8CsVfiSOUlqAzuGQ8K08lP3nk5sM0oF1tlwqYcDk87y2fHflwivHVbnyRgLGnuYN4tFKD7+o96ZQ7oBesFgs5sN+sjiOyCqb7/rywwDMbXu5+ltM+CtD2XaDgZ1Jmgz5wVaSHVCax3ElXagXwCNSZdhcbqhpv87SGZpXA5XsKIw+irM4gSEkdsdWuOSNbBezdIoShWQfpiGYFf3Kv0cRfoV67t5Ohz0NSIaYx20bhO7tAEFGBIOXQOdW5OQBEsLjsKlKITz61JO8+9N3Ynpl5Xv2KB9yl1AkVovwZbDwquKX8K8/Bfv/C1Ce+g1z34Yv/xJ86b+p9DXgl295RJFdIkN7cYwOinRlEvSJHFuZI2mEgyplGkHQdFtREeNz3h7aUyZX7ejhr5x3YkqHdyVuJ500+Z3pD/KLxZvJJAz62pPMF2wMf5p+YMePk8QhtXAMXvO7MHQ5ABPH9rONGQbn95G9/J2MyiE6ypMIK8//SHyDd5p3LcmkUZ0gYulEFH6tLB1QC9muNw6q9FoIbIlLPfXSGxQL2K5km3OSS/OPww2/oAZnOcxJz4sOykZ7JeH7udVOohMPg3+Rb+QlxgFuKt9DojzPPnlhYCdVw/U85mQnnW4W5kfgP3+xMiBcD4fvUES/60a44odUTCSShLClfJJRYwfFzJYgaNthqzHitKv1AqVEJwWZxsue5vvFc4z0XMfW/rDKaCLTVeHhjxtDzGZ28RLjOQCmUzvZ6k3QP7cPgHRWvUBJdYR7ByRrK3zPz4YC6KaA4c9mtaXTbU1wWg5wRG7nInE6iBNUlEOxq+5T11ZmLbU/hoi0m/xUoPCnAoVvwCOf5tdzH8OLJAUULVfZg5e9Gdm1jbeaD1cSfiGyGVFuIkhndnNhzI354xi+5bPHOsSCbCfbdRFc/y648kfPLcIXQrxJCPGCEOKwEOID63WeqKWTThjh4oYtVwRZOjo6/pSnsgb2CpW6pRXf8MI+6NoOZhr2fxm5oBSK2z6MedkbgnPlk8rSmTX6g9/laK8k/M6tmNKhl5wijrahwLO/ydiPRRJ6dkHnFsSEyox5IaEIz5s9xpBQUftJ2Yvo3QOoDB5oXC2zIg+/tAAHvw0PfQpQRLerdBDaB5ntvZrsmCJse/YkCenAha8G1JS5I52gV+RICI8Oey44ZMIU2I6qkri1cAiHBIfkTjrTCbZ0ZbC69nBMDnOxMUbaFOxwRxnyJv0aR6bKciiq65je8Vp+oPwJ7v+xh+GmX4d+9cK7tnOWN5lqIdTM7jczIfvotKbpshWJv0icqChT4HhSkVs0Dz+nCX+ipocPqjbRK41ncIX6uzGrCP8Kn/CHhLICf9S5FUck4Lqfh1SnsnacMpSyFEW7UvhRS0e3K6m2QLyl/CqmZA8fKH4EgH3eRXX3kHV8S6fTy6og8TP/rlJWo8iehr+5tnIF6dRBGL5K5RX27FIvnUh7tliK8MvpAQZRaZndziwuBrJN9eOhzgxjsh8xsZ8LxRjj3d/HNn//XtMQGOkudZ3lLC4GRZlirP0KdYJML6c7rmCbnKA/q14AmUU/Sy3VEcwSnVQ3WIvB+hONqE3ZKyJefH4KSlnSbp5x2c8RuZ2LjZDwo7EcYVUVIOscZtHyyMs0hh0h/NxkQPjz2UU6KZBKmJCbJIlDm5cNPXzLVUH4vr04w9dymRiptHSK4dggNxl4+DIfIfzRx9VqbB9jcgDHE/DDH1OW3rlC+EIIE/gE8GbgSuCnhRBXrse5ipEyuJmkGea6br0SrByeY9GF8mhH5BCukWLYUDdZBWEkw9mn4cJXwaVvgP1fCbItvK5hRPc2nvEuUOdKKktnRvQF58/RVqm8O5VqGhILyqbJDELXNgAuM0Y5wbCS453DwXTvcEoRfltuhC2a8OnD7N2BK0VA+LUUfsJaAGSlh6+XbY88rDqzKdhlHYFt1/DoYh/WnDrvFtsPbl+m8p4vFqfoSCfoRw2Q9lJYBVRl6SgFNFw4yKnUXmwSwQYYV+/sYVz2s1XM0C8WSeLQJXNkkkYQLE0W1fGSPduYope87Q+GrmHKIsOFxgRvMR9hqv1ish17mJB99Dgz9PgprReLU9jRInWeVATs2WCmlIL0FVW/yNEuKxdFRRX+K419TPa+mCnZTWL+KLgOl3EcUCm4tuvyJnk/h/pepZ5pSpE4hRnwbIpGByXRBlYep1iZWqcJP0sHP2X9PnP04JhtvCB31q2N73qSOdlJAidcxDfyUOWHTj0Oc8fgLr8ei+eq5IFBZckE2WJ6fwYrz5B9mtHEbkX4IovjSbrcORbNPjIppYKHutJMyD6SIw8AMNl7Le2pBD1tSTpSJiLdEXj4RdGO48Fou+qzDF7KYtsOdjBFuy8QzHl1H0l1hKWkE/4irKpZS5Tw+4j8LT8dJF+MyQEOeLsYFnMk8+PB/Qrgz7xkh7p+2bmVbMmmQAbDLlQofG3pvGPm0/x76k/UC8m/X932TLhQ055RL8+e3VgDV3ChOE2GSKZbsVLhpxIGKdNARpS/HH0UgNN+DGJC9lWKw3OF8IEbgcNSyqNSSgv4IvC29ThRqSoPf1jMIoUR2ASiPE+Xr/AXZTtW25aAVBOm4AIxoTzNXTfCFW+D3DgdJ+5QB/eJ+h5UkKqYGsA0BNMRha8snUjn61SbPGwTM/SziNU2CF3DKpUTOOwOqxo7XVuDr4yn9pCTGToKp9iCGjSTspdMJsMYAxFLp4osclP8/ANv5A8Tt1Q+UP+FBRJe+BYZw2WHfQKGr2Zc9tNjT4PnMez4hH/xD+KKBJcYp+gyHdqFmu62lcMMiZRpYDhlXi2eYDh/gJMpNVvS3uU1O3oYkwMMeTMMoWZZ3eSCKqZlxyVZVB88jsEAACAASURBVOon06vua1BeQQjGzW28wrqP68VB9ve+mrLjMUkf3d4cg7aacSWER082zLJyPBn6wj071fXOHaeUVC/kjvwoUWiFVlqc4UrjBLNbb+Ko3E5q/ihMv0CbsDiaupSUcEksjjIoFpjq8nVK2id8PyGgaHRQFG1gLTI/G6o6F0Ol2/k4Infws8aHuPdln8MhUUH4dz0/yZs+di+2q9JA5/BJ8eSD/s+HK9ofZDIdvp1v3PodFUdyyyHh6z6V8xeFnX4SA49Dycuw0oO0izJeKUePO0vW7AuE0lBXmnH6EdLFlibzvS8CYFtPRgUjU10q7bE4T9HowPUkp7TCH7yUfNv2ynbOaUunMyjPYfkvwWqC0yKmrz0ZlnEAP+VR9c/Tsp8H5NUA7JxV98bxJB0USeAo2waY7lAzxVNuD4sltQDOdGpbOtfY+7hYnCJlEtiAPU5Yb6jP8u9h727KA1dgCslQ6VjYvipLB9RYMP0XgUxkEDNqJn27ez0AY7K/kivS3cHakfXGehP+DiC6hnzU/92aQkpJyXHpl/NQnCedNNjGLG77FuhQC6TM0jzdvoe/SDtW2xBbNeEbBjea6qGw6yUqMIdg4MS3cKSB0amO8WVexz87r2Oy/RIMIZgmVPh5MkssHYDLxAiGkJTTg2AmWTDUd47JbcwV7ODFAFBI9XNSbmWgHCr8KdlLJmky4m1ht5gkZRrcVP4ePPIP4blOPkBCWrwrcSvmw38X/n7ev/Vt/XDgG+yVp0hiw/A1jMl+9e/CDNvd0+RlmnL7MPNtu7lYnKJLhiqrrRQSfpfM8TelD/JPqb8i5eR5rkOll2mFf9XOHsbpp9ebZaurvteHCtqm/OJaqeIkWdlOZ5citmgBtVPGdga8ae4T13FP/zsp2S4Tsg8DyUX2oeBzA7kXgn+rBTd+e3v8mkZOibFuRQ7pXGUZA51WJ3zS9vov4pg3TCZ7FO/UkwDBdbVNKJVd7PC7bapD/fS/WzI7FeGXc4jSAie9ISSCIpklsYPTbi8TnYogo5bO/YeneX58kXzZUVk60idFnR0y9pTKBtGYOQLpHgqiHfPBv8Gb9O/FkL/6u1rhjyh77HDqcqyML1Lyk/R5cywmBwKhNNSVZtxXoc/JPRgp9cLa2ddGd1syvPbFMYpGO44nGclcol5Qu26k0L5T3XpMbBLhiynZHu4OlvCvrSpwqwl/uKeNXnzCb+tThL8QKvzR5AVMyj72zKlZj+e6fCf1Ad6X+LKybYCZDlXmZEb0kSs75GkjYecrLZ2SQzslLpCjJIVLlxsG+rvdcP/dQcef3fbuptinxGN/LuyHgaWTaAvud0c6EZTTyPeEi87u9JRgnKC/cr+DTA94zobUxF9vwq+VMlLxGhNCvFcI8ZgQ4rGpqakVnaRsO/yccSu/9MSPwld+hUzCZKeYwu4YVp0GMK2FCoVvt21hyF9UZRqClxj7KZrdakbQMQA7rsd0S0zSSzqpyGzGHOL3nfdAIoVpwDRqcNhGBodEFeErS+dKQ6XDldID/jHUatyjchtzBStQY54UlJO9HJbb2WafZIuYp0SKLO0kDYNjbOMicZretgTvtL8G93woPNfJh3GMNN92vx/jzj+Gab9DLoyoeMT3/TQcvZvrympqyfA1nHL9l9XiaXbIcU7IYWYLNnPpnewUU7RHVnqmC+MwfRju+CN+7divcIk8xm9av8Knb7qbfT2vUpfrK/yrd/QwJvsxkOwpKyLqEQUypiTtWzrp0jSTspcef8P5aDnZW1I/yc1bPsDvpX+PrJeh7HhMyF4ALneeZ1L2kpVtbIkQfqXC3xX8frTjKgCSC2EOOkTWFfiKLtWzhWNyG6nSDDz2T0zJHka6VAZK55Qi/FKnf1xdF2ZRzTZKZgcFlKVjlLNM0EehfTt50V5Rghvw1yQooo8GbUfmCsHftaUTYNcPqFxuP/MHUEQ6dBn3dbye18pHWTjiF36rtnR0XaLRxxhL7KCY6MVuGwyuvdebJ58cYGdfOwlDcOFgh8oeAp7wLglq0v/uW67gQ++4JpzdLI5RNtpxPUmJDD+c+DRc93OUOhThHzX2smD2h2QYsXQsTfil2oS/rSdDn5+izOBlgaUjEUzQR0c6yUPiGi7IPgKey67yQXYZU1wgxgLCH+28Ck8KThvDLJZs8mQw3ULYR/JT5C2HK8VxTOEvciycVBlsqM1atMIfcjXh76LQsZuiTNG3GCX8WUhkVPFEX+F3phMky7OAYLZdEf6s0c8T3iXkM1t50rtoqaUDG2LrrDfhjwK7Iv+/E/z13T6klH8vpbxBSnnD0NDQik7iPX4Lf5z8HK6RhGP30mbYXGMcJTdwdUD4ifJCkKWzSBt2+xaG0Apf8grxNEd6fiBMc7nk9QBMyP5AASWDXF1VQ2aaHkBgJZTyqciPT3dhGRmuFD7hp3zCN3zC97Yxm7cChT9HJ4lEkkPeDrZ4k+wRE0zKXkBVpDwmdtAncuzJ5NkrR5Xy8RUJIw8x1nklf+C8G5Fshzv+SP1+YRR6dsCNvwjA2+ZvoUQKBi5i1Cd8uXCKCxjlqBxmJmeRTQ4yLOZIWaEfnSqMwx1/CPf/La6R4le8/8VXvJfjJRVRQKjwBzvTDG5XsY5dhXA5ezd5tUGG7ZEpTzMle+lMJ5ZUzNzv7WLfwJtJpxKU/FW6Ez4J7ZUjnJKDHJB7GC6Gg871vJBAencHv59I7mROdiH2fRE+ciV86hXw0KfCVZV+Fctk9zBHpbKXjNOP8ynnhyi3qb7YP6tW01pd/nGrFH7Z6KAg2sAtkyjNkJUdzHddyhzdFSW4df/QdYWi9YVGZosVvwssHQjWanAy4uPPHoP+C3kgdRNpYdP2zC2K5NvUi5FMj3rR59QqT0Yf4VDyckxDqFgSYOQn6WeBQrKPGy/o54n//Xp297cz4Sv8x71LgxfjRUOdfN+u3jB+kR2jaHRiu2o9hGNkQAicru040uBJ9wLyydDuJNUZHKts1lb45Qjh94ocEqFKMuSnYPwZil17cEjQmUnwqPlitU7h9FNcU1Iipp9FDCsPwuDZ7lfxWuuvOCmHyZUc8jJD0ol4+FaOcjHPNUZozfRmD4T/dmcDQt4iJ5FtfWo8S8ELche9i5G9DwpzagbduVUp/Nv+gF+2biZjz0NbH9O+wBuTg+Ro57Y3fpe7vRdXWjrnEOE/ClwihLhACJECfgr42lqfJHvZT/Dr1q/x6It+H6xFdp34LzpFibmhG1WNG8D0PXxPCvJkcDu2KhLConP+AINigUNdLwkPevEPAirAotM2dUDUFAJDCLUUvWOQsqFIoGIVqxAsJvr92h9Q8BX+lE/4x+SwKqDmq7EZ2U0yYXBY7sBAcqPxvE/46rwnhVJPLzefDV5cTD6nCnCNPc2prmuYF73w8vfB899QBLEwohRv/4Xw/e8hic1RYy8eRkD43uTz7BGTHPD2MJO3mE8M0idyJPygWE5mSBYm1EKWq3+Cz7zoFr5rXxW0S9+T6IYpv/2OVwOwPR8Ooi6ZU5aO45IpTzNJLwlT0JGq3BKuaHm0pfxNbGyXku0FhG8gGWeQ/d5edpSPYKAX3xAO5t5QX8yJHkbFVnWfeveoKfPdfx6kmaZKynM3u7YGhO+1D/Iv7g/i+IQ/kHuerGzD8IVDtYdvJTopoJbJpwvjZGnnkSt+lz8w37fE0oFwr1/L9XBcVbBtZFY9T23zzEcV/p6XqUqsT39RpZvaRZU23H8h+4wrmZOdZMozoboHlanT5RPQ/AnIT/F84nISpgiuKzF7mCQOhZTqj92ZJEnT4D7vKg5e/B7u9F68dEWpftnZecpmu1+1NBwXyVSa99q/xcfKb6Ps93eMJCRSgVgqmf7LrJHCJ4ed7Faz5Pw0nHyI/NbvB5R6fjL1YhULe/Y/ubb8GAD9YhFh5yDVRa7sclxuY2qxTN5yKZAh6RYqAsVGYYqrjaOUpJplds+HfbXHm8P1JELATqaQ/qzRcjwOeLvpWng+9NuLs2qldOcWda8f/QyvKX2XNmcO2gcY81S/Oe6oF6AWRs65SPhSSgf4NeBW4ADwJSll/SpGK0TJM/m6dxOLW1Sn2LlfpSFOD9wQKnxrgS6K5GhDYgS5x1vEHL2nVB2Xg53fHx50+4vJZbZxUO4IVuLqaakmOldKlVli+oRflS65mOgnIXxFl1QD4Hupl/PV9A8zR5eqa9OlCb+HlKkIH6BX5APCT5oGJw1F+K+0HwhPMPmcyuTwHEY6rlEpmS/5VWU7PP0FX+H7BPjK/0XB6ORZcSmW6zFDDw4mHFJF4Z6Te5jNl5n1FYk5pQbAC3IXqcmn1XR15w0kTBFup+jveAWVxdb0eoPovq6d3gLXLdzBy52HaS9PMyV7SBgGHelExUbmetPxTNKk5G+mMUMPrr9d4LQxyAtyF2lZYrtfX8jxvIiHvzM41gy9fDL1Lvjxm+Fd34KX/SaU5ukpqFlXuzWDTYJERx8n5VbybdtZfOnvUCKN0daHJU1M6TIqt5DR5aIjPjYoxZqXbepP9gJZ2c5CcgvH2LFkwRdQsVlMyfFYKNos+jMcfR/m8QnfTKkX1Zv/QmXh/PPbYdInpv4LKbki8IUrCB+UkFgcV0XxgAOJyzENA7dN9cP0jEqdLKYHg68kTYM8bTxy8fsokllaJCwVvojKZodP+F5A+OmEyXe96xhjQKUhA/hxAF3fvuzPhgOFv+9L8G//jZ13/Hd2iwmu2tHDy7YLzM4BFX+TLhRnKW2/UV1WOkEp2cdD3W+Ehz/JJfYLeAj6xSKmnYd0ZxCQPT6jXqQ5mSEVEL5qa6o4wzXiKPd7yv7pnFW0VBYZ+rxZFQxOJdghpnG71RiyXckBuZtUeS58DsU5xTGdW1Vw2S7Q7c2z1z4C7QOMOIrMR6S6zzpPv9LS8Wdmm53wAaSU35JSXiqlvEhK+WfrcQ5djVF274DePaTzpznsbWcx0R+8PZNWlm5RIIvqgJ6furWFebpH7+Z5LmDejExDDZMvv+Tf+Vvn7YHCDywdITCFqpvNNT/Js71qSXo14WcTanAVZJqyqUjhGeNKvrz1NwDBbM6CDvXimaGblGlwXA7jSPVYNOEbQsUP8jLNVUU1gGWyAyb2B1P9kc6rlBuVaocLXqmIfHE8VLwdA3zsss/zcfEzWK6Hh8E0fRijKgPkgLebmZzFrKHaLPxKkC94uzD04Nx5A8nIyi4jqvAj9fNp61NBLMJUtA53kVeP/QN/IT5O0isyKZXCb0+Z5H0/W0pJwXJ8wlf1kEq2ausU6l7MJLZw3FM22AUikppXHbQFpmQPL6SvgqverlTvLkUavTPKD2935vh/7b15mCRXfSV67o0t91q7uqqrWt2tlrqlbu1q7QtYCC3YIAsEboExZjAYY5AxeADDjC2/9/yejRnm2ZgxxsZvGH+WAAPCLLbM2DIw2A9hSQi07xKqXquru2vNLSLu/HGXuBEZuVZmVnVVnO/rr6qzMiNuRsT93XPPbztBBmCb3Mn41Wv+HvN7uYSSdiwh2QE/ZWOB8ZYavmT4Vg5LguEDPASzIqJtHK2XriQNejewYsVTcg4QMHwPBhaQ4ZKGYfLd5m2f4+GYUq4bPh0Vz8e3ReRHVfZjkJASw/S/A1YWz5NtMCkBNW0cZkPIvcxj+8upwODLSBrpX6gpAywXO/CyEa7P4DGoRd/RJCw/Iw0+XyTke4pUi9J56K95EuDBh1GY/h4+ZX0KKeJhZ64CIzuiAi4AwJ28jH8tx4RtUNw19GtAYQoGfDxMz8EQFmBUFwE7qwz+S7OccCwhDccXBl/kwgyVXsZOeggP+WdiFgVVUuSAvQODPvc9ZG2KSXIMlRz/TNXz8U3vCp5L8A8f4ix/+bgw+GNipPx7TvmHgMwIXijx50UWoJNlPtarht8XyAc0ZVFg25UAgPv9s7kuaJiAU4Bd5Rr+AuMGnwnt/HR6CJkjD+IH9MKanpiLSMOFyZsxI2ApBiWglJdRxZXvw30jtwOorXEzLxaQGTagNDvPZ8qgLZRdwLThDu3Ec2yCN9eGiZcYX4yOskFYBgEhBKZh4Hk2AYtVMMMGwCYv5gz/8a8B4+diiRaCsgpnXCdC2ViI8ZZSY1jwLbV9PowREObhBMvhMIZxfKmCY8Lg48hjWGaOGgvMFLD5nBDroyTYzmf0ejWEKJb/pM+176w3h3z5CLIi1HOGDcKkJNQSruL58BmfFCnTQMn1VKldKeuctMbwPOP3TibOcaetmCyFLbwvgZnGSdcJyyojZwKpARSEwS94JzBHB5WmX/aCbM+0RTEjFtyX2abgOBGGXzWyWBQMHwDmWUaVXjYNonoMF4SDWi8dXap6ymEr/y9xmGwCxs8Lxr7nFu7AlWWkh3eg4vq437wY/819HZ4efhVCyI9zDX/6h8DkRagwqnamv1V9t8oQL0cYPhA40WuqQjqBb6FsBgxfttXU30/zY6HrJY9dhcX9C898G/jGHVjY+ko89ob78MTlf4jz6fPY9uifCCM6DGTF2LJjMEZk+K8J26SY89PAm/47vmG/Bg9YF8MiHlLFI4AdMPyfCqlsCQ5sv8hlJFGb/8oSv44PsTNxmA2B+nwhPuDsxLDPI2w2OxVkSRmVNJ8DFdfHcRRwcN+HeLeqR76sSTpinuy9FR7Es5IdwQ+XJ/C37rW4T+zEpPQZaseYGPz2IJ1dKcvgBbgA/MA/GyeWKvjoPY/AdwZgVcMMH3luNF5D7wdhHh40zq9pYi6zOaWWKVmKQflEVhn9wihFGf6cwY3UMQyoMCzX92EYBBk7kDKO7r8Xf+reqiaFlHVmMKgMqqnJPU/7U3BHzuLFoA4/Aux7B3zGgn62O68LBqExXsugqLp+YPBlCJ6/DQDB7GIFMyLyCAuHMEdyKmoDWy4EDCvU7tCgRMlcIUkHCAw+4wY/v/AcTFbFMcZT9WfJEAghyNqmaoJSqkhjG0g68t4eFeOYtzdjBoNYRrqW4ZspXrguNQDkNqHoRtobUgpM7kP2KI+8GcEc5o0hxUwrnq/uYcY2cYzJ7fimgL1aab6gVJcBKwvDtLDIUuoU88iqaBuTEkUSCqLc9ILO8Kue0u+BcHjqb9CPcilHxxXv5T9Tg0BmGBXXx/nbN+Pj7n78r8MRf0FuM5cbDj8CTO2D5/swhcH/vn8uHrnxi/iM+1os5wInt2T0chx1NXwAVSOnSmyYmqQjYQ6Mhz4jnxvXF6XEX74fKEziN7zfxO9+6xkcnHg1vuZdic2P/3+crGSGA4Z/2uVICUKRT3GDX3F9YPJi/HHq3Zg3hT6+NA04OUUgpCN4iaVBwXhUlijVfJn/EJaRxgP+bkUm4AzguD2OPJbhoIKtFt/ZFlN88ZJkYGHPW4BNZ/N690XhtJXz7LxfwLEMzwNg6REcWPDxn8l7MM34d1GSjp48mRLlK/pQE399GHxX62y091YsXXoH/sm/CN97egZ33f9TLNI8nIpk+IKNZYbhwsA19BEwauNJ86xwXQ6IbvYGVYbU1CUdqeEDsWneADBnCEctG1R/kxMkbRlqYnlWDlWYiiE9Iwz7UTakJBTTIHheGnw2hcrIWQAYn/zn/QJ8nwVlFYZP51UdgRDD59Uyg+YRB4VD6QlhlGeXKpj30yiCG7B5FHBELgCiqqO+zZfXAYgz+Hysz/iT8BlB/gTXjH+/+hb809gv42F6lvicoRjlcpX/TNsGHIuiWNEZPmfbi844AIIDdFwxfGXwJQNNDQLZsVCPBIWtl8I5/hRyWBalroeUYau4vnK8p2wDM0yXdMRxCAlknVQBlkGwyDRJh2VQEeWjKSWwxPWRDF/X8IuVMMMvaob2gD+s/E8KZ/0sv69Cr694Pk4bTmPvlgL+5+NHwu+VjNN3galL4coOZeIZOZ7fhT9wb4dlWeojknBIImJFGb6m4VfMrOp8Jp87R/NZpAaFwbey6jsBwjnrCAN3w/+FmbKJhZKLsuvj8+6NoG6RZzGnhzlpoBZw+itViGvOMVV4L8DvfdHk9ylVmQXsfGgXBQDL4nmGV+GLiDMAEz6eyFyMKkwVmYTcJiwKGXYTOYktBmfcy86mYOwAbMsC9t7K5TLf5fdp+9XA2+8Fdt2I4wWepDdPC6h4PvZMBPWIlKSjk0vT4RJowvBbQ0lKOqYBpArwrvsdFJHCU0e4rlu1CrDdeW7wBcO3DBMn6BBM4qM8fiE8I1XL8F0vpEsGTltu7KQEJHcCUUnnpGT4rBC0LRQTL2MbaoLLhUZOiid8Xjtnmo3CEAbWohTTJjfez7AplIZFSvvFbwPsDDwWKauw81W8wp9m8GVpAznOA9Lgi/MdXyqj4jMcFwvVAs3jOTbJ66CfeYO4bsH14AxfavgRwyoY/kE2ijlkkZ7lPoHH2TZ8bfBt8CmfhBk7kHTk9cjIKB3XUwz/AX83nvBPgysSh36KibDBL80HBn/sbGD8XJF9HRnX1D4QMJxPn8Mo5rBsDoOK71Fxg2uTsQyl4b+sa/hAwHSdAkxKAxIBYB4ZVDWGL41mISUlnSjDL6r7VpK9YC0a224T1AB+8avAz/8ZANmL1cCr92zGQz89gRm9m5k0+AAwdQlcjymGz681P1eofaVk+GWv5m+h7w3ANbNcwxfSFRDW8LMj4US1gOEz7ps443pgzy1YqrhYqriouD5+xM5AZUj4IjJDnOX/+v3Axb+MrG1gJGtj+2gWttai0vV9FK3BYIxOLpTXwcegRT05eTCxc/jpMJd/FcPPjmFRBFeM4SQmKNfyFy0uLcndn2VQVYaEj3WYE4FtVwCEYG6YZygfrvLvfs7kQDAUk8ISPSVC+MDjwKt+F73GujD4E4Np3HbxFEZFo2754E2f4A6xxexWjC4/iyGyqDR8wyA4KUojVLdeFWjy4Abk2GIZZdcPsRarLsOXkk6E4YvjH8OAusHSEGTsgOHLuS0n2D/4l+JN1TvxHJtUi4xpEDxunYOZkUvwPf88LI2eB9z8R8DVHxDHjRRO+5mPcuNgBuzTjmzZn2NbwEDwMNuJkayN40sVvHBsCQsWnxALNI85Ywjkt6e5IxiBHwPgUTp1Gb5wFh/ACE6wHAyReXiAjWKp7CoDkHVM5bTVe9CmTBGlI177mn81bq78AVION5wvsnFsJTMw4QoNfyFgjvvvAn7uv6qInxAm+U7lWvoTOMRF0eH3SMoEclKnLAMP+rvwLNmGl9lYSK5QoZmpAkxKMO/rDD/L49MZE1FMwuBrGr68TSUh6cjmMEVhqBzTCGdi6hjZCYxyPbvserBNilfv2QzGgPue1Fi+LK8wtB3IbRINa2hg8Ku1Rj1g+HWcttRQ5SKqZk6VEVYMX7tG+RFRZkEafFlLx/OB2+8Gbv8iQAiWyx6WK55YaAlK53B/GERBN4zsBKgB06D4t9++Dm+4aBKOFRh83wdKdjjmf7Hkhu47TWl5DU4enmhFOidq/R+WGfO5TYHBJyexSTQlWhCvyWfDNikwfi6vTqqPVWBh7FJ4jOAZj9+DvVsChm8bVPWUCCEzDBgWeo11YfAv2DqIT7zxfIwVOGu0DQrd9k1P3AjHW0aeFLEAzsYsSnBSMFl361UwKVGT7G8feBnXfvxfMLtYDuuSWlgmlVE6CBh+tBPVjDmGCky84I9HGD5FWmP4qv2iWKh8UEznzhPnDDT8ijOCH1z7PzDNNqHqA7jsXSrZhjGGECHLjvJCcBqUU04w6u/4F+CfX30vnmOTmBhM4cXZZUyfKCI1wh/kRcojh/SLaWm7CEK4IcvYRq1hPe8X8MQr/xzTbAxzIsywYhWwiAyWKl6g/dsGliMMX4/SiTZsl3HMz3ubYRGP9wCWYZmS4RMCEIJixasNjUwPwh3ZjRsoj9+WCXFy9yMnomNRfJddiDfgE6jACu8UdIZvECx4epROBmVXNIihVBlNpeGXXZVhXKx4mD5ZxI7RrPi/XGxojTwYhe8zVD2evbxnooDJwXRY1pEMf4pHJrmahg8EBl93tJpRp22U4Wvf3bXycH3e2zgapZO1DTjZQe6cVU5bUVrBY3zhMEx1rqWyGzD2c/cDkxcDU5cgCsc0QAgRrTMDhl+2NOnLyWGp7Km2lQBgpXWGX0BpaDf+3d+F9MhWWAYJJB2d4ZMTGMVxLDEHS8JmyDFack5Ilp8JG3xv9CxcXP4MflDlWv7eLQHDtwyqekqsBtaFwY+CEBLaXh4cvgQLQptTDJ8SHDc3ocQs+FOXqNZmAPD0kUUsVzw8dnA+LOkYgdPW0HYE9Zy282QAb079Gb7hX1mj4WdsU2nWspa9zrbGRVlaU0k6fFegMhYjhlDXUutBfjbQOAkOGTzhaGKAP9SUAJsneabsojEQctLy8YQlnbdevg1ff+9VgcNYws5ieceNAICTjE/6UoazvuWKqwxAxjGxXPXg+0wZISnpuD7DYtkNtVuU0UCSPW0nh4PEKydgUkC4gqoOf3IfdlBuHKspYfCNMMPnE5OqxdEJSTrCgDh5mAbFom9BhuPNs4z6HgYN7p9k+AAwKH4/ulBGxfUxNZRW4wW44Qw1Z49BRWObhBBcfvoIHj2gJTPlNvOonr23AuDPh2EEBr8U45iV9yQIy4wz+Py7e5Zg+F4QLCCv0VDW5gbxqt/g0UUInht958LDcPmirnYchTHgnfcBE1qEUnQIZsDwPZ8BVkYlULlmVvg2gnlupLTnwslj+sr/A2+u/CfkHAuFlBVIOrkxlKwhFJmN08hRDHnHcYQNKVmx4oWJGS54M5enIiGxOcfESeTx/MwSCAF2bc6rxdagRPWUWA2sS4MPIMTIqj7Fo0M8bE1q+Cal+ObA7XhL5aMw7QxMg6h2aUfmeaGq6RPFMAOijZ22UQ3fZwwn7c3wwUsKA3xLa1CCtCbpyOPo51IGX5xzMGNjU94JIkqiBp+1YPDNMMMHeE14ANgizrdv+zDSI1yOWTYKNZPeikTpZB0TZ4zlgo12ggAAIABJREFUEQc5VlkqoCzimZfKXiDp2AYY48YuCK81FDOfK1ZVA3QAyoA/Jwz+DnKYJ/Ecfx4oTITOX6x6KrpDB9kasMeq2N4Hko5cfInY9Qknri7p2IGkY1EClwWvLSCjvodBqXK6Sw0f4PcSAKaFw1YutrqkA2jN2WMgnzV5jQtpM6xdUwN4x7dVVy1XEA35jJTiJB0qGb74W9Rpq313z+Jx+GGGz8et7td1H1MlSuR7dGZbdoNGIyeWK6FjNIKtOW1dn8G2DByXz5jB57dk+PmUCU+rWgonj6WKjypMZB0DhbSFQ2SMO17HzwUxDDzDJrGLTKPgHsNRDKnFSDlt5TWbvAh434O89paGwQy/1z984ThGsg5sk2IkZwfRfqKnxGpg3Rp8nZm7vo8HB2/iBZVEAoRpECzZm/Ag2w3DEBKNuAeH54PKhHoBrFAcPiGQi7Q0+FGG7zOmxlHD8C3NaSs+ZmvGdHNBMnz++f/n9efi47edpxaFqMHnEkJjg+9EtuxAkAg0MciNzo17x1U56GVjoGZbH47Db3w+abTnhaTjKoPv1oRzLlVcNbFkaQWAL0hDusEXBnwWBcwhi8vpE9j80jd54avz9qv3+aI1YRzDN7Zepn5naR5yZ5sUZU3S4VqrXOAjerbU8J0CTEPIL04OZcKL6OkMX16vQjrYpUiDIJOutgzye60YvrhudfsPQzM+SkYxsVzxwBjDk4fn8bcPhCuEeh536suFVka26SRDGngVpRPV8LXv7lp5MMZZrxGRdOSCpkPF4WvfSQ9DPSFaDcaeMwLbMAKG7zHYBsVxEe5bIty4nzbCd5U5x4QfcdouCqd0PmWikDJRNbPAh18Edt8MgxA8zbZiN30ZuSov8icXx2D313iMe7cU8Ee3nYezJvK45kxOKDblHfU5kwYEsN8wm7/l1ESI4XsM06kzcYv5aTxa5g+GXgfGFFEakuEfntMMvj4hpKRDCAwaMPOyeiDCq7bns6C1m67hG1GnbYykUwgz/E15rhO/JNLFo7sJLuk0viaWyd+g166RDP/qM0Zx497NuOWCLUCxAoDguDMJy41IOjTM8BtBsrVFIoyEcHItaTq2TERZLnthDV989mSxiq3DAUMLErwIvkB/Dr+KL6L0xAu8ebfWEFwatDiDT8fOwgJLI4sSjwZBfUkH4M8S0Rc3qeGnBmD6hH/GzmGZ8mtZqgYMX0k6qVpJZ/okv5fyXherYdYeF6lzcrnCDW2EbWYcA57PUHZ93HX/T/GVB6fxxn1BDoZk+DIsM4jS0fwzRvj5qKvhE4PnI4A/+4ZYwOS4hzO1zkeZu6L7JvSd5onlqpKnmkGXdFyfwbEojjNRapvwcW3K2cjYBvIpC8wKfCxw8lgUkmbWMVFIW6FFz6AET/lTuM34HvziIo6wc+FWA0JHSVjWjAMhBG/ctzV0/UdzDg7PldX4m/loeoUNwfBlY4nj1gQs4SwyKQkMuMic9XwG32dK0okexxCslFKttAIaM3wZuhiN0knbZk1YZiMNXyIUz6zB0xOv6iDqtAWCuPDJwTT+/K37MJpzeF31Dz6JlzLn1Eo6ejZlk7kpr92SwZ1WTNQkWa56atJIJ+xi2Q1r+MKwzxWrIadwWivh8EXnDXjOn+DFw/a9PeRc1uWhGlCKn+BMHEceqRS/zo5loKzF4VtmwPBrjqFr+FL2sbNYEmUDpOZraMZhQNPwcyleJTRg+FLSkRq+yMaMMQq/+cWH8YEvPRzL8AHOmhdFXLsO1/djo3TsGEknYPh1DL6Th6E9h3IRMQ1+/DiGL/9e9X0cmS9hdrEcYvgnlytqB9oMMg6fMR4WahtUSTrLIgku51gYydnIO2YonBROQT3/WdtUReMkCAGeZvw5pczFUTYYknRir0kLOH00p0ibSUnD3Vsvsa4Zvm1QVZXQE/HvKYu/xhm+kGhEEbCq52N2qQLXZxjJ2phdqoQ0RRmhYlAo4+r7TDG6GiMsWLdpEOWEk1E6jknU2KJhmUDA+gwafsCiko50HPs+C0or1IE8/qJm8OcEw3ei0Sz5cezafLImcSlUS6fJ+eS1WzIHgSrAhrYBOAnGAjapG6pwWGbAch2TIusYKFY9ZDTjS0wHH6q+C3869c+YOO9NoXMreSjO4AP4LH0ThssHcJXsg2xQVFxP7Zwsg6jr5US1bFuTdCiPmWf5zTh+nIbObRi0JvEK4P6AlEkxV6yCEGCswA1BSX3/WgenxOOH5rEp74SctkCw81kqu1go81BV1/NhGhTLFRdVj2EgbQVO25goHVkbKS6CR2FoBzD8U7UYll0/tNP7wKt34eozRms/Bz5/XI/hfXf/CJtyDt5xzQ71N8nwW4EqheH6cH0ftklxQjD8JaQAMGQdA3smChjNOTBMA2VmwSFVUXqBh1vmUyYuP3047JcinOFLHKfDGJHz2/NbHmMU//HG3eqamyIBcjWwbg2+Y1JsH83g6SOLqHpMbWkzQuskwshTEjzopWrA7q88YxTf+PHBkCEMZdoKY+cx1oDh8/dalKqsRAAqDh/gbFclXpnBgycZvhWh0crgi3P9n996HI8emMNw1m4qscjPxhn8uO37R24+q+a1aGmFRpDX7n7ncuB1nwPGzwfwXX4cFb8vDFXFRbHigRB+73RWnbIMtfDoNXsc08CDbDe+c8ltuN0JO46V8Yxx2gLA4+bZmCmejleLv9umNIyBVGLUY/haHL5qzv2zf4w//usHgbmAqRsk0Mx1hu9YPCx3qeJhKGPDMY2QoZULZVTSWSq7ODJfRs4xayQd6QtZrnih0gKmQTG7yPXxEe0ZiXPaAvy+lN3gOa3Bdf8J8D8C84eH1Tn05+DXf+aM2s/IYxsUrudjZqGMquerBC+AM/xWjaleCsNn/DtISWeepQAUkU+Z+PO3conv4/c+iSWk+HU1bfX8Zx0Tb71iO956xXZ1bIMSHMYw5lkGBbKMOXMUmUpA6GJlrhaQtoNn2DZIEqXTbbz9qh1433VnKubO2Q43tGHtPnDEej5T+v01gqWEwzKFl50G5RakXMR/D09Q3w8cZa4fvE9G6QDcOPhqIQjONZS1Q0ZHQj5wZdfHQqmKL/zwZbxwbBmej6b6p/xsSNIpVUM1X5ohVFqh2QIjndxWGjj3NtU5DECN01Zq+Bmhl+tG1jGp2gnoOw5HMeFatiRllXoMX45N/j0ahy/jpeX5wx/W4/DFGNJjqrqpYnI02CXIOHxAMHxx3hHhkLYMoiQOyfCrke/1oqj+KMtGy3EDGsMXce1AIDXOCofoSM6ukXSiBl9eF57LEnN/DQuw0pqk4zVd+CUsg6DqMyyVXcwtV0PBA50wfLmwmpTgOOHXfk6EXeuJgLZJscRSYIIULJVdkfFaez7+TBM8xUQuijmirlXV61zS0cEZfmLwu4rXnDuB156/RckpUkpJWUEsu6nFJRuEh1keEgz/ip0joCRSEEqTdOTndB0yehNlqKRBqaqgKI+jGH7FC+rLU71UgYl8yqx5wKTxqXo+/u7hgyhWeVcoP5p4FYM4SWe+WK01aC0cA0BTCYmKmGNpwMJJPiIOX5Mi9No3esJUmOGHDScAeDGTp5mkI8cSsK5oWGYjDT+opaPCDf3ACMtzU83g51MRhi+OOZpz1PmiMosXIRAvHuNO3uVKoNHL51NfOBeVwefHO77EnYXDWVtz2sbLNirvo0kkSj1JpxFk8b6lsou5YjXUB0FKd61ALkpy7hkGwd+zK/GNnXfiKOFRV7mIwV9ECkzct8Wyy7X9GMiv8rTPdfxFZ1Rdq6rHOpZ0dHASmkg6PYElVlPpLHVMqjF8qhlxri8emSuBEt555/3X78LF24IsPslKdUmn2MDg+0xOeu60lUzUNCjSltyCu0rSIQRqR2BQgnzKrHHaWopZ+fjKQ9PqGF5LGr6MwtAZvhuSG5rBbEPDB/hkkwYzbreU08MyK556b4jhW4Zi+JkWGb5ygMY0IQHqMHzXV4aUx+HX0fBPfyVwyTuBsT2wXuRlkl2tKF1RY/gy4UbWUKl6LLSAjYhyILZBa522kW2/ZPjFqlfjtA0zfFF9VOxyjglJZzTn1Gr4dcJuY2PwNRgdGHyeYepjqeKh5PoqPFKiXYYvFwyTEhRpBo+M3IS05pBV7zcolpECs7lzfLHs1pYCkd9LPNNf8l6B6y88AziYjzhtW/uuzcavz8F+YkMYfNeTDJ9LKdLInzs5gAMneaSEQQl8xnB4voRNeQemQXHHq86MHCsc1QMgtA2POm25IxUiqSue4euSjnQeS2O3d3IAm/Op0DHlw/7s0UU8emAemwsOjsyXUap6bUTpBBOtHWalXwMgaP/bCLoeH4r5ppLhh5228rqkIs1D5Ou6pCONVbTKKRA4tet9NxmimrHDBn+hVEXOMUUPgjoMP7cJ+NlPAAgMn+sFZaflcCTDlw5o26Coel7omkiGz30IUadt+Hu9cEyXdMLfL3B+65KOZPjc4A9nbUVKjotEp1wqbAKUwW+yXQyF57aw8APc4T8vQiI9n2FGi4YD6oSBxkAuiHqCm4yaWyy7yNpGaC7YJsV3vPNx9o6tyABYLIWzt3XIz/2YnYGjl/0yMt98XC0s3GnbPDGsGRKG30NIDV/WA09bhprIb7pkK950Cd+6GSK87sh8CeMD6fhj6U3MxfMkdch8yqpZtWVkEE+0YIqxGTQs6UimLJ3Hkn18+s0X1YxBlVA+ugiA1xH6x8eOYKniNs1SlMZhscwLeEnD1CtJRx5bGjDpJPe1JDHJfJWkY8VLOpk4p60V00xCIFTZMAbSuKQiGv5csap2PIGkU//6KKetz2pCIU2xW9PPsVTx4JhBmKkM1bMMimI1nG0aDct88VjQMlLmTyiGL5zfi2UPi+KZlDWeZhfLSFl80Vws82MeOFFE1jYwFImZl9+nmfHVWX10F1oPpkFwUiw0AHBQ+Mvkzqd9hh/spKQPrlT1ati7bVB8yns9br/8OmTA/Vb1drX6M21QgpRtqMAG7rRdOcNPNPweQtaAd0VGYM6p1cUBqNo4h+dKGC84MUcKQiT1KpGSZeQdMzYOX0YDydBQACIOPzD4MoGLirjtTJ3IEvlZQoLJv3MTjxhZKnsth2UuVdzQlreVdHZ1/jaidABequH8Ke5Q4zWOjNBYgKBE8rIm6TgRp21GTOJ05HWgVusGUBO2GIU8f0bT8Muuj7nlqsqENSOLQhyCKpCsZodHCcFY3lERV3IsKU3D15220V1JnKQjb3E0ukrez2MLZbWQ607bkawDoj23PgO2DmdqHLOmYvitPUvye7YCk1I1bgA4eLKIjG0oA90qe44afKoZ/IWyW7NriYYyzxWrocxnHTSykPGMeL6Adstpy2vpJAy/J7BEhIwnMvJ+5ZrT8bPnTdS8TzL82cUSrtw5EnMkPQ6/VtLJpUxUPR5rLyeRL0vkihssGZshwkMBoFh1lc4sa/TU0xcBqGqBRxfKIATYLiotLpVdbK6zUKnxa7V0BtIWKi5ntfV07thj6BO9BYP/J7dfGPq/bVIUq15o4ciJEsmlqqcMYA3Dt2IkHWUY60s69Zhq1GnrCEnnZFEz+JLhNzBE8ntUfR9lr5bhf/CG3XhvNbz4OGag4QeSTjgMNfq95ktVHFus4PTRLJ4/tlTD8OUCclSriS91+tnFivIV6MZZz2CWaFXSCTH8NqJ0jswHBv/QXAkZmzc0OblcbVnSke8rVgMNX9a2ipNroqHM80U3lPmsQ78+0TLm1RXE4euQvozVwLpn+LLWSdVnMCnFni0FvOrszTXvMwjBzEIZC2UX20ayMUcK2A8lQREqqe/Jh0zX5gJJJ8LwjbCkwxTDJ6rPayPIh24s76hwv6Wy25Rx6yVqbYMqJtnIoNUeo31mp0OeU3f+8knFnbZ61Iw8vGNSnDaSwWjO4R3IxOuWmOi+yLjU49blfag3QZWkY4YlnZPLFQymuXGU97vRgqjXiIkyfEP4YwZU6YFgcVFhmcppG1zLuNIKckd3tqitPhcx+FQYp5mFQBeXDP/4UkUVNNON89ahOIMvo3Ra1/BbWfjlMU9qDP/AySKyTiDXtUo8ogw/qG3FQz5rDH4kO72hpKMNwaQUKdsIJVZ2JSyTJqUVegaTEtFyzm/IREyDqAly9kQh9j260zYapRMYfL38K4+80UND+edpKA5ffkTW96npHhWBNAgTA2llOJY0X0A9hNLoTapkk3YYfkjS6cDgy8mqSwYZx8SiaIQhvw8hRBnjlGXgLZdtw3f/4yuFRCakNeGsc32G2//iB/j9bz2hjllxG5T4Fa+nrCCfwjYoGOMRLYUaDb+xxAaEY8Il6uVQOJYRG5YpIc+nP0tPH+E+m3NEbfWTEYMPcGlMZ/jSaTu7WMZI1qkZ09bhWl9Vu1E6QOsMXwYvSFRcHxnbVHJdO6UVgPA1lyUu4iJwLI3hVz0fyxUvlPmsQ59DhpB05MJS8Vh3JB0zKa3QM/BCRb7S8OtBv9F76hh8FbNPESvpAOFJKtsOWjQIDQWgqmXKz+thmTfsGceuzVp1v7jvJB66ycF0KC69qcHXJrFtULgmP29bTls9LLODZ1+eS78XOYc3QSlpUToAl3WKVd7ERJe6DEoAjy8aUrt9aXZJVVwEWmD4Jg1dO/m+uVhJp/4XlQugjHzKpUycXK7WfEf9HI5Jkbb573EG345h+E8cmkfK4jtUOU4gvIhnHSNUB6pc5fVmZpcCSYeQwHF+WqykI5LNmjL89qS96HdUYxb9i4HWwzKdOIYvfHDcPxVeoB2N4UsprJCqE5apfRdL+NqKVb4Lr0RannaK2I5XfcK6N/hy5dd7b9Z7H8CN6EBMtT8gLOnIZ1dKOjKRQ69iKWPjpWanR+mYBuXxwRFJ53deu6fpd5ITY8tgKuTEbC7phI2KdBb30mkbhR3jtM07Fg6cmK9pWJK2DJxAtWZ8pkGAKjS5jKFY8XB0oYwlwfD0mjhx2JR3VAlqPq5gPLKapcq0bcDwo8lsOae5wU9ZFK87fxL5lBVIWNr5A6dtYPAfPziP3eMFtZOU53AiDF+29QR4xdAl0WBkRCsxbVAC32ONNXyz8b3tiOHHvC8jmpIDncThS4ZPlYZfrPihAnv6+yuur8JCW2L4wuAzxuWxqse6EofviBDcitsdn0A7WP+SjsiglD0960GylLMn4pt5ABGnLanH8INJyiUdzkKjmbYARJtDNyTptAL5kEwMpBVTBJpXr9QnnG1SJZm0w1o6ib/WEWj4wWcv3TGMF2eXsVwJG/wgCSs8PtX2kRIYonFNqcrbCj52kHd9aua0/a0bduNvfiWoi68vKgHDbx6lIw3wrMhm1fXj6P10FMM3sHs8j3e/Yqf6W5ykI3VexhgePzSPPRMFdX3mi9Wa8gdZ2wjtCspVH8dlHZ1c4NCX45KdtnS0HIffgbRnhr6jjC4yFCNv2WmrJB1uvKXEKsMyo9nVIYMvGH49DT8UpUODaKpixeua0/aqM0dRdn3c9+TRFR+rXax7g29p+nlDDV/8rZ6cA0QZvtBuq1LD5w+Q3teWLzIi+StSSweAigDwtbDM1r6TZPjpkDFqJunICB+ATy65lW5HwyckKCvdSu3yKJTB1yb3DXsDJ7rOzpSPIcLwDU3Dlz4ayeh/Mn0SAN9pWQapO8a0bYQ6aekTeUA5bUlozHGQpYCPzjc3+IGGX3s8vXBe4LTl3+nQXAlzxSr2TOTVjmCuWFt7JhPRrsuuj2NiIQoxfEIwmrNDkpaEnAdtJV61qO3p7HiLyHXRNfyWGX60tIIWlrlccWvCmvUoHSmF1WP4+uIVCq4Q/ZW7YfCvOWMUY3kHX37w5RUfq11sAINPQ+WR6yFg+A0Mvsbwa5y2DTR8KTsEDD8ICdSrZbZqQOVDNzmYbkvS0T9rmzTEONuBXnCuXcQ5baeGMjhnkl/3dCgck4Z+BucXDF90KtNT9B85wEvfVtusbBiSdKIafgOGL5nizKIw+Jo2HGW+doPrbccwfLlbfFzsWvZsCRh+nMGPatelqqcx/MDgU0pi5RwgcHA2u3a6ht/qZZaLSNoyVBezrKMx/HYlnWo48apY5XWpakp6a5FUssNbvbBM/btwSUeET1dczBeroZpIncI0KG69aBL/8tQMZjQnez+w7g2+KQqXuc2idFow+GN5XotkKGPXOG11Df/IfAmM8WYqMkU/GocPcIZf1Bh+q1tjORknBlOhh7sV55meSRlX46YV6J2/2oU0dtHF4sY94wDCkzWQnGI0fAQa/mIpCPV7ZJob/IrnN4000aEbuNrEq/rHkc4/yfCzjRi+Wf94OqOOhmU+fogb/N3jgcF3ReMPHZKxE8LvUdn1ldSk72ZMSmJDMoFAtmymVRsdMPygQqqh/CRZx1TjbpV4RKtlSolVZrpHF2i9wqxs+FNP0tFJl8zMB3j0luszNe6V4raLpuD5DPc+eqgrx2sV697gW9Jh6jV22u7dMoB924ZiIxckrtg5gv//I9dhfCCljF00Dv/l40Vc+Qf34fvPHgsSryiticMHgIxl8uJpYlPQaly7bVLeGDlrh2LoW/m8NCxhht+uwRfS1gqidKKSwc3nTsAyiOr+BLTA8IWGLx2mO2RSUqnK9dY2GL4TknTaSbyiyDmmin/PN5F0CIlnz/riFBRP48/LE4fmsX0kg5xjxiaeScjeAlnbRMoyUHa9oDRyNtDwP3jDbrztyu2x36cjDb9lhk/EOE0VGJG1DTXu9iWdcOKV3OlFNXxVP9/VGH6dTNtoaQUp6ciy6e0UGmyEM8ZyMCjBIXHcf3zsMJ45stCVYzfCujf4KstVJF7Vw2vOncCXf+3KhiyZEIIx1YkqYPiWVi/l5ePLvDDUQjlIvBK1QvQoHUA6bcNhma0gbRmYHEzzEDtKglDHFj4fNvi1ZQxagc6w20Wc0xbgE+DfP3Z9qFtSXOVM/byyTtGCiLzYK0IWp48XUW4zSSYcpSOTlJo7bQFuBOTWXNfwo99RLrJx0p3qHUCDom2yLeYLx5ZwxhgP1ZWhqNExAwHDzzoGHJO3bJwT5a/1heIXL98WqgKrQ5VWaCMOv2WGr3U5k4YzYwcMv1WDTwg3xHNFzWlLCRbL3JjrgQz6cSvieujMvd734iVMgnl9qMsGnxBeDXeh5IIxhvfe9RC++qMDXTl2I6z7sEzJ8Jtp+O1COt6OLZbhmIZiL5JRVUU3HkJ47Hq0WibAJZ2DJzVJp8Xx3fGqM5WRA/jCUXb9liQdW9No9UzWdqAYfieSjiykFmOMo71QVdGxyHtVPwMx0WUBOylblF0PVa+9KqC6b0MVe1NhmY2PM5C28MxRzs50SSd6P84az6vEqXrnlxVTgYDhn1iuhOoRpS0Di2W35rpILTzrmChXfZSqHgh4K79WYWuSXyOEo7VaO7ZcQHOOqRbVrGMEvo02FujxgRSmj/P+ALKc+JJi+PVLK8yL4nj1/GWUhsmMZPgyv6FbBh/gfoSFUhXFKn9e6/kVuon1b/BFkkMzDb9dnDs5AIMSPH1kESNZWzGiWeG8q3hMSTpGTMcrQDhttQYorRrQcybDRiNjGTiJakuauq0xfHm+fhp8WzPWzZCyuHQVNZyKhRk803ZRLH5SXy1VfVRcrz2GL947qBmDViQdgBsB6WDVjWv0O0bb6ekIZXGLz3miD/KJpSoGs1pPXGnw60Tp5BwTBLxJiuezuqWA49Bq8TR9wTZavM7ymBnHwEBa9jcwQ4ttq5gayuD5GV5uQmZcy+egntNWxuHXi9ABgkg53cEMAIfmeH5Do8+2i3zKxHzJbepX6CY2gKQTRMh0k+HnUxbOFYbXMakyGCdE+deqG+wqrJo4fP7ejMjik693ookDQd/WljR8Mwi76zxKZwWSjphArZTUnRhIq2bu9c7Pt/LC4GfCDL8dAyLfO6gl3bXitAXCE1U3ru0siPquRZ5XlgGoeD6Gtd2PlCzqRenkHJNLOlUfC6Xa6pGtjKMn9fBlZU/HVPcq6xjqmrVzvyY1X49BCCgNkh6jco2sMFsVYZn1smz17xJl+IeFU74XDD8IFe09/173Bl92vGoWh98JZFVNR2ubGJZ0eJSOIxxoMmQzeJi401bPtO0E8gFvLUoncNQGpYjbewykQeqolo40KC2sbr/6itPxjfdeXfO6rrNK/wgQGOuy67fdnUgZ/HQ4mgVoruHri4Qu6bRaJx4I7othBJKO5zPVvGRIN/h1oqskw886JhyLoux6DZt9xI8j3j8QRWe1dISkY5sqLDPnWMqINisaqENPGpO+HImowZf5JzLxqiHD154tINgtHBYMv14WfieQGr50JCcMvwuwDJ6Yw1jrzqVWcZXW6FxOFDlBXZ+JRh/876Wqr8Iy1cNkGShVgz6qKzX4reysuxGlI3XelUTptGIMHdOInWBqwdHkDyAwvKWqJxKv2pd0dGPQjoYv0SgOvxF0Z7Zy2vpM7RiHsjrDN0OfkQgzfMoZftlVSYGtoBOG33ItHUl0HANX7hzB/33rubh42xD2bingM794Ma6qU5Y8DrrBNw0SGkPUaQvw570sonQaGXx5z+Q9kHNrZqEMSvhi1S3kUxbmi1XMLTfODegm1r3BNw2qenu2w7hawcXbhkTGasDwZaKLKuegefpVr9NIA28ZXtYJYwYCFtKWhq9n2rYr6WgGt13YyuB3/ugFYZk09J0DScdvu06JEyPpZCwDlCA2I1WHbkAahWU2gqVdU72pyglhDIY1DT9dp9iYHqUjwzIXy9W2nLaW9nw0wooYvmhC9ObLToMhomFuOme8rWciyvD1qR2tpQOIfgee37AWPhCQGHkPZM0rn/GFvdXFrRUU0v1n+D0TjQghdwJ4J4AZ8dJHGWN/36vz1YPVwYPZKlKWgWt3jcIyqJp8C6qfKF9k9LBJmRgSbY4tFwLSqYavlRRuhnA0SmeJV/I6rqwefuf3Qg9G0X2pAAAYlUlEQVTLDDF8MWHKVa/tOPxA0gkm3RsunsKu8XxTSaQuw+/A4Js0qPfv+b6qADoYI+nUROk4QZSOIxjtYslt0+CT0M966KxaZhCHv1JMDgb5MmbkOYgLubQ0SaeRYaURDR8Q/ao9v+sGOZ+ysFhx1aJ+Sht8gf/KGPtEj8/REKEerF02+ADwp2++CJQQLGjZnkBQi5xqDF8a/CCOmr8udyArl3Ra0fBjnLZtaviq6UZHYZnCabuCexF2cAbHGVoBw09ZBgxKQkXG8ikLV+4cbfApDn1X0KiWTiMEjUyCekVVTdIJO23jE5Ukw8/Z3GlbqnpYjGkI0gg9rYdPZRx+ezvKOIzlHdULN7rwxxl826SYL1ZR8fyGztEgAiw4Xsbm0TTdNsiFlAnGeKtHoL3w2U6x7sMy9W1itxk+ELDr6ASRRlx2PQKgMgHlOKIMv2NJpy2Dr2v48aULmmElUTp2ixpxIyiGbxDll9E7hfFStu1p+CnLwOfffqmKvGoH0hBQgrZrG0nIhVivU+T5DCeWKiAkLBvpTdF1SIafS3GGP19yUfVYW1E6pkYIGr5P1/BbjdIxgyidlYJSnpX90uxyjcGP27HaBlUBFY0kHSMSkgsEC2w3QzL1cUyfWEbOMVckc7aKXp/hvYSQnxBC/ooQEpvaRwh5FyHkAULIAzMzM3FvWRH0rWmr8cKdILq91hl+VNIxIpEQJVETpEN7rx7IVj6vx+Hv2VLArs25unVV6iGoGtreOAE98arzxVdVdBTVMgFuaHl1TOG07aCy4dVnjnYUhSENvmMaoYW/LYYf8YvINngnlqsYTFuxDNY2wgv1eCGFj9x8Fl5z7gRSlqF2B/kOGH4zOYxSou5/yz1tafcMPhCEZurlyvUuZjpsk+LYYvPQyiDxqjbqp/uSDr8O0yeKDUNFu4kVWUBCyD8RQh6N+XcLgD8DsBPABQAOAfgvccdgjH2WMbaPMbZv06ZNKxlOLPSHsRcMXyLKiJSGTwJGtqjV/gBqGX7Hkk47TlstMmfX5jy+/ZuvaNvIWQaf7J2UR7Y1rbpT6E5jNdFtQ7VFLLt+19rRtQIZymlr+RhAezs2XaYCRP6I7+P4ciUUkgkEzv7ogkYIwbtfsRObCyk4Jm/ZCKDNOPzWGD4fa3vOe720QjcgHbeyAQpQ38FumxQvzS6DEGDnpvod5eTzpBNFOb+6LumkJcMvdn33UA8ruvKMsetbeR8h5C8AfHMl5+oUnTKudmHQoHUcwB2H8nWZuLNco+FHJJ0Ox9eRpLMCY2gZtOPFqaBirjvXcXWjqDN8gO8gOMPvTju6ViANgW3Sjn1GqrSCFgfuCklHD8kENIbf4PvpfpnOwjKbj122mmzd4EuGv3INH+DZtvy4gaRTr0aO/F4fe83Zqk1kHOTt079TpkcGXzL8uWIVu8frN17qJnoZpTPBGJO1P28F8GivztUIeoJPLxk+wB8qyex1R2zgtA1asgGBdi7LvHY6PJV41Ua1zJWwX5PGb5tbwaXbh/H5/9CZVi4RitKJxEunRIZpt9rRtYJ8ygQhfNek68nt7IBqGL6osHpiuYrJwXC2scysbrSg6X6Ztkor0NactvpYWzX4Z47lsHU4jW0j2ZbH0whv3DeFgbSFQspSu6l6WdE/f8Ekrto5indcvaPhMWmcht8jSUf3JfQjQgfordP244SQCwAwAC8C+NUenqsudK24104RO2TwAyMua7HIEgDyWVIaftUD6VAiATorrbCSzj2WQTp2MFNK8IpdK5Pu9PLIKhvWDhh+2eWJV/3qF0opQd4JasJYBlElr1uF7HilL2aS4Z8TYaT1Mm116H9rq3iafD5amCtysW3V4J89UcD/+tB1LY+lGSYG0qrMs16fKg5vvuy0lo6p+1AkeiXp6PelH0lXQA8NPmPsrb06djvodZSODsukgGhgIw2/LuksVVxVdhUIa/idSiRAe5m2utO2U1gG7ak81gx6eeZgKx/4JpYrvD5RvzR8gMfJ60ltVbAmnwjDivg2MraB40sVnFjuTNLRy0G0Y/Cl/NPKZ1bS+azboE0knZaPE8m0BQJJZ7CLZRUAhLpn9Yvh929GrBJsPUqn55JOcHwVpUOJ2l4vld3Y8LFS1e9YzgHaq6XTDYO/fTSLbSPtRfZ0EyoLUsu0VZKOZajS0f1i+ACfsI5miNt91vQ4fICX7fi3Z2dRdv0ap62Kw2+woOkMvx1J55LtQ/ibX7msJcltJeG53UZQ+2ZlHDY+8ao3YZl6Ke5+FE4DNoDBN/us4UuENXz++mLZDY3B7hbDF7VDWpFZWu1Z2gjvuHoHvnXHNR1/fqUINUAxwlt5x6RYEI0wVvId28Xu8Tx2iEWwkx1QlOHfsHezqv44FGGWbTtt22D4hBBcdcZoS/Kiug8reHa7BUrCO71OESfp9MppCwQsfz1o+GsCZh8ZvjQwaVHHBOCTQTK/UtUPxdtK5l+qrMzgpzpw2vYrgqUXCJKDAklH7451VHSf6ifD/8Qbz1e/y0Jd7SAah3/p9mEMZSycWK7WSjotNP2Wz5Zt0LYT61rFSjqfdRvNonRaPw7/2Q+nLcCzbWcWyn3T8E/dWd8irJCG39uvK881nLUDhq9p+EDYp6Cctq7XR0mn9TjrtQpdSpD3VHdkykYYq/UdO2H4qqicCjmluP7szQDCDciB+rV0dMhnrh123y6MNSnprMzgy52NoRHF00ayGEhbGNXKbnQL/Wb4p+6sbxG9rqUTOpeIcOAGP4jS0dsJxmn4VY+tqAqfTDZpJQpxUz4F26Rdd0D1E4au4Ufj8E1Nw19Ng9/mjs2KMHwAuO3iKQxmLGyPhDGePVHAz1+wBRfV6UsLBAy/Hf2+XejRUqsNGtnpdQp53/Sii689bwL3f/RVKz52HKRz/JRIvDoVoEs6vY7L1hm+HqVDRHmFUtWP1fCBzrNsgfYSr27Ysxnf/9DP1PSPPZUQ7XgFaBq+qGwItBZL3gt04rRVtXS0Z/Sy00fwo//86ho9PeuY+H/3X9jweJJM9LIgl9xddbNkcKeQhjqzQoYf7FqCZ0dvZt5tSEOfMPwuQU+86jnDp5zJ6+w5qPFRa5Rl6zX+vs7PO5i1kLYMbMo333JSSjAW0zbwVII0inpZ3JTG8CVWi+HbBunYaRtt0tNpbkZfGL6x9hj+isMy+7xrKSiG3x/uve4ZvpRZgD5o+CZR3YYklEEyDQDV0IOkM/+VLEaFlIV/+8h1fWMJqw1VjyemtILuL7HN1TFElkHbNhjKadulIcsonV4yfPnMrmR32i202o6yGaIdr3qNQhKl012Y/WT4BkUhZYX8BvKU0hBFx2CLjlydMjmJaCTHesbrLtiCwYyNfMqqlXRCDL832/BmsM32S09QkTXcrTacqX5q+H0yjo3QLNO2Vcjb1i+Gf8XOEbw4u7TinUmrWPcG3wqVVujtTRzJOlgccCMGP8w8orsMxzKAkrsiSWejYSyfwm0XTwFADcPXd1f9qqUTRdYxO5rAnewM6sHpZ5TOGmD4NJKA1/FxYjT8XuKVu8fwyt1jfTkXsAEMvtnHKJ3ffd0eVF0fn/nuczXnVJ2ejFqGD6yNSXMqIhqdoW/p+xmHr+PDN56F5arb9ucsg4TCAVcCufC1UymzXayl0gpymq/YaRtTHnk9Yd0b/BDD7/GDKfW4OIZfr5erZGIrlXQ2KqLx12GGvzoG/7QOy068/qIpXLZjuCtjSFkGbINiNNc7qW8txeFLRr5SDZ+uoe/UC6x/g99HDV+dUzf4EQYap+H3c2zrDUYk8Uqf8KdaNvGdr9vbtWOlLANffc+VOH1Td0oRx8FaS5m2YggrL57Gf66FyKNeYN0b/HAcfn8MgC4lqDrdiuFHNHzVELwvQ1t3aKzhn1oGv9s4ZwU9B1rB2mL43XHaBk3M1+ezsz6/lYZ+ZtpKhBs885/1GL6MKlkLoW2nIoKJLmoEhcIy1/3jvapYSxp+t+LwLYNXsFyvIc7rnuGHa+mspqQTxI7rCMrirv6kORUhr58sL5HSwjI3OsPvNdZSlI7sk7vSEgWWQfHN912j+uWuN6x7g2+IbFbG+qjhxyReSSZfy/ATSWcluOmccZiUYMsgn6AJw+8f1lI9/JvOGcddhcuwuQtZ5GeM1W9yfqpjQ8wIi8br572C3nQlmngV3WUohr8GWNKpiELKwusvmlL/D4VlJgy/pzAoASVrI8IsZRm4cufoag9jzWNDzIh+1+1ulHhVn+Gv/qRZD9CdtgnD7y1Mg/aNRCXoDjbE3Qq6CfXf4EeLe0UnSLS1XYKVIXCCrw2pYT3DpCR5bk8xbIjbZRl869kvx6gVknTCiVc1cfgJw+8qpHSWsPvewzRIKM8lwdrHunfaApxV93Pr2biWThKW2UvI65lE6PQet196Gs6b6m2sf4LuYmMY/A7qk68EcZJOc4bfp8Gtc0iGf6pl2Z6K2LU5j12b86s9jARtYEMYfLuLVQhbQXx55PjiafUWggSdIWH4CRqhWq1ienoapVJptYfSEVKpFKampmBZneUbbAiDb3axCmEr0BtvNK2lYybF07oJuYAmGn6COExPTyOfz2P79u2n3JxjjGF2dhbT09PYsWNHR8fYELOCa/j9u7mhpitKw69XSyeIKkmwclBKYBs0YfgJYlEqlTAyMnLKGXuAk8KRkZEV7U42xKywOmgqvaLzNQjLrMfwE0mne3BMmiRdJaiLU9HYS6x07BtiVliU9DVKR5d05P2pWw8/CcvsOhzLCJW3SJBgreGee+4BIQRPPvlkX8+7IWaFaZC+9t1sh+E7iYbfdXCGn1zPBGsXd999N66++mp84Qtf6Ot5N4TBt4xVlHRUPfz4OPygxWGfBrcB4Fg0cdomWLNYXFzEv/7rv+Jzn/ucMvj33HMPrr/+ejDGcOjQIezatQuHDx/u+rk3RJRON5tDt3o+CRJ12ka0ZVndMZF0uoesbYbKJCdIEIff+8ZjePzgfFePuWdLAb/72sady772ta/hpptuwq5duzA8PIyHHnoIt956K77yla/g05/+NO6991783u/9HsbHx7s6NmCDGPxtI5m+RsHopRWiTcxrWxzy1xNJp3u483V7kLY2xKOd4BTE3Xffjfe///0AgP379+Puu+/GRRddhE996lM455xzcPnll+P222/vybk3xKxotuJ2G3GSTs4xcen2YezdUgi9VzL8JKike7h4W3cagSdY3+i3XQCA2dlZ3HfffXj00UdBCIHneSCE4OMf/zgOHDgASimOHDkC3/dBexBokpiZHiAk6YhfDUrwpXdfgVfuHgu9V2r4iaSTIMH6x5e//GX80i/9El566SW8+OKLePnll7Fjxw58//vfx9vf/nbcddddOPvss/HJT36yJ+dPDH4PEJJ0mhjyRMNPkGDj4O6778att94aeu0Nb3gDrr32WlxzzTW45ppr8MlPfhJ/+Zd/iSeeeKLr598Qkk6/QQiBZRBUPdY0Okgx/CTxKkGCdY/vfOc7Na/dcccduOOOO9T/8/l8z+LzV8TwCSFvJIQ8RgjxCSH7In/7bULIs4SQpwghN65smKceZKJXM+KeVMtMkCBBv7BSSedRAK8H8D39RULIHgD7AewFcBOA/0YI2VBxclLWaSrpyObmiaSTIEGCHmNFBp8x9gRj7KmYP90C4AuMsTJj7AUAzwK4dCXnOtXQao2cpFpmggQJ+oVeOW0nAbys/X9avFYDQsi7CCEPEEIemJmZ6dFw+g8ZqdPMkDuJpJMgQYI+oanTlhDyTwDiUr4+xhj7u3ofi3mNxb2RMfZZAJ8FgH379sW+51REq+UcVGmFxOInSJCgx2hq8Blj13dw3GkAW7X/TwE42MFxTllYBmlJl5f12xNJJ0GCBL1GrySdrwPYTwhxCCE7AJwJ4Ic9OteahGVQtJooZ5s0kXQSJEjQc6w0LPNWQsg0gCsAfIsQ8o8AwBh7DMCXADwO4F4Av84Y81Y62FMJlkFbTqbaPZ7HjtFsj0eUIEGCjY4VJV4xxu4BcE+dv/0+gN9fyfFPZbQq6QDAV37tyh6PJkGCBGsFH/7wh7Ft2za85z3vAQDceeedyOfz+OAHP9jzcyeZtj0Cl3QSnSZBgjWLf/gIcPiR7h5z/Fzg5j9o+Jb9+/fj/e9/vzL4X/rSl3Dvvfd2dxx1kBj8HiHR5RMkSBCHCy+8EEePHsXBgwcxMzODoaEhnHbaaX05d2LwewSTkiTUMkGCtYwmTLyXuO222/DlL38Zhw8fxv79+/t23sTg9wjtOG0TJEiwsbB//368853vxLFjx/Dd7363b+dNyiP3CJaZGPwECRLEY+/evVhYWMDk5CQmJib6dt6E4fcIdp8bpydIkODUwiOPdNlh3AISg98j3H7pabhi58hqDyNBggQJFBKD3yNcumMYl+5IeqsmSJBg7SDR8BMkSJBggyAx+AkSJNhQYOzULcq70rEnBj9BggQbBqlUCrOzs6ek0WeMYXZ2FqlUquNjJBp+ggQJNgympqYwPT2NU7XZUiqVwtTUVMefTwx+ggQJNgwsy8KOHTtWexirhkTSSZAgQYINgsTgJ0iQIMEGQWLwEyRIkGCDgKwlbzUhZAbASx1+fBTAsS4Op5tYq2NLxtUe1uq4gLU7tmRc7aHTcW1jjG1q9qY1ZfBXAkLIA4yxfas9jjis1bEl42oPa3VcwNodWzKu9tDrcSWSToIECRJsECQGP0GCBAk2CNaTwf/sag+gAdbq2JJxtYe1Oi5g7Y4tGVd76Om41o2GnyBBggQJGmM9MfwECRIkSNAA68LgE0JuIoQ8RQh5lhDykVUcx1ZCyL8QQp4ghDxGCPkN8fqdhJADhJCHxb/XrMLYXiSEPCLO/4B4bZgQ8j8JIc+In0OrMK7d2nV5mBAyTwh5/2pcM0LIXxFCjhJCHtVei71GhONPxDP3E0LIRX0e1x8RQp4U576HEDIoXt9OCClq1+0zfR5X3ftGCPltcb2eIoTc2KtxNRjbF7VxvUgIeVi83s9rVs9G9Oc5Y4yd0v8AGACeA3A6ABvAjwHsWaWxTAC4SPyeB/A0gD0A7gTwW6t8nV4EMBp57eMAPiJ+/wiAP1wD9/IwgG2rcc0AXAvgIgCPNrtGAF4D4B8AEACXA7i/z+O6AYApfv9DbVzb9fetwvWKvW9iHvwYgANgh5izRj/HFvn7fwHwO6twzerZiL48Z+uB4V8K4FnG2POMsQqALwC4ZTUGwhg7xBh7SPy+AOAJAJOrMZYWcQuAz4vfPw/g51dxLADwKgDPMcY6Tb5bERhj3wNwPPJyvWt0C4D/wTh+AGCQENKTbtRx42KMfZsx5or//gBA5yUUuziuBrgFwBcYY2XG2AsAngWfu30fGyGEAHgTgLt7df56aGAj+vKcrQeDPwngZe3/01gDRpYQsh3AhQDuFy+9V2zJ/mo1pBMADMC3CSEPEkLeJV7bzBg7BPAHEcDYKoxLx36EJ+FqXzOg/jVaS8/dfwBngRI7CCE/IoR8lxByzSqMJ+6+raXrdQ2AI4yxZ7TX+n7NIjaiL8/ZejD4JOa1VQ09IoTkAHwFwPsZY/MA/gzATgAXADgEvp3sN65ijF0E4GYAv04IuXYVxlAXhBAbwOsA/K14aS1cs0ZYE88dIeRjAFwAfyNeOgTgNMbYhQA+AOAuQkihj0Oqd9/WxPUSuB1hYtH3axZjI+q+Nea1jq/bejD40wC2av+fAnBwlcYCQogFfiP/hjH2VQBgjB1hjHmMMR/AX6CHW9l6YIwdFD+PArhHjOGI3B6Kn0f7PS4NNwN4iDF2BFgb10yg3jVa9eeOEPI2AD8H4C1MCL5CMpkVvz8IrpXv6teYGty3Vb9eAEAIMQG8HsAX5Wv9vmZxNgJ9es7Wg8H/dwBnEkJ2CJa4H8DXV2MgQhv8HIAnGGOf1F7XNbdbATwa/WyPx5UlhOTl7+AOv0fBr9PbxNveBuDv+jmuCEKsa7WvmYZ61+jrAH5JRFFcDmBObsn7AULITQA+DOB1jLFl7fVNhBBD/H46gDMBPN/HcdW7b18HsJ8Q4hBCdohx/bBf49JwPYAnGWPT8oV+XrN6NgL9es764Znu9T9wT/bT4Cvzx1ZxHFeDb7d+AuBh8e81AP4awCPi9a8DmOjzuE4Hj5D4MYDH5DUCMALgnwE8I34Or9J1ywCYBTCgvdb3awa+4BwCUAVnVu+od43At9qfFs/cIwD29Xlcz4Jru/I5+4x47xvEPf4xgIcAvLbP46p73wB8TFyvpwDc3O97KV7/7wDeHXlvP69ZPRvRl+csybRNkCBBgg2C9SDpJEiQIEGCFpAY/AQJEiTYIEgMfoIECRJsECQGP0GCBAk2CBKDnyBBggQbBInBT5AgQYINgsTgJ0iQIMEGQWLwEyRIkGCD4H8DrtODqO3P/Q4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(range(m), A.dot(x_true),range(m),v)\n",
    "plt.legend(('Ax','v'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Recovery\n",
    "\n",
    "We solve the relaxed maximum likelihood problem with CVXPY and then round the result to get a Boolean solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "ECOS 2.0.4 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS\n",
      "\n",
      "It     pcost       dcost      gap   pres   dres    k/t    mu     step   sigma     IR    |   BT\n",
      " 0  +7.343e+03  -3.862e+03  +5e+04  5e-01  5e-04  1e+00  1e+01    ---    ---    1  1  - |  -  - \n",
      " 1  +4.814e+02  -9.580e+02  +8e+03  1e-01  6e-05  2e-01  2e+00  0.8500  1e-02   1  2  2 |  0  0\n",
      " 2  -2.079e+02  -1.428e+03  +6e+03  1e-01  4e-05  8e-01  2e+00  0.7544  7e-01   2  2  2 |  0  0\n",
      " 3  -1.321e+02  -1.030e+03  +5e+03  8e-02  3e-05  7e-01  1e+00  0.3122  2e-01   2  2  2 |  0  0\n",
      " 4  -2.074e+02  -8.580e+02  +4e+03  6e-02  2e-05  6e-01  9e-01  0.7839  7e-01   2  2  2 |  0  0\n",
      " 5  -1.121e+02  -6.072e+02  +3e+03  5e-02  1e-05  5e-01  7e-01  0.3859  4e-01   2  3  3 |  0  0\n",
      " 6  -4.898e+01  -4.060e+02  +2e+03  3e-02  8e-06  3e-01  5e-01  0.5780  5e-01   2  2  2 |  0  0\n",
      " 7  +7.778e+01  -5.711e+01  +8e+02  1e-02  3e-06  1e-01  2e-01  0.9890  4e-01   2  3  2 |  0  0\n",
      " 8  +1.307e+02  +6.143e+01  +4e+02  6e-03  1e-06  6e-02  1e-01  0.5528  1e-01   3  3  3 |  0  0\n",
      " 9  +1.607e+02  +1.286e+02  +2e+02  3e-03  4e-07  3e-02  5e-02  0.8303  3e-01   3  3  3 |  0  0\n",
      "10  +1.741e+02  +1.557e+02  +1e+02  2e-03  2e-07  2e-02  3e-02  0.6242  3e-01   3  3  3 |  0  0\n",
      "11  +1.834e+02  +1.749e+02  +5e+01  8e-04  9e-08  8e-03  1e-02  0.8043  3e-01   3  3  3 |  0  0\n",
      "12  +1.888e+02  +1.861e+02  +2e+01  3e-04  3e-08  2e-03  4e-03  0.9175  3e-01   3  3  2 |  0  0\n",
      "13  +1.909e+02  +1.902e+02  +4e+00  7e-05  7e-09  6e-04  1e-03  0.8198  1e-01   3  3  3 |  0  0\n",
      "14  +1.914e+02  +1.912e+02  +1e+00  2e-05  2e-09  2e-04  3e-04  0.8581  2e-01   3  2  3 |  0  0\n",
      "15  +1.916e+02  +1.916e+02  +1e-01  2e-06  3e-10  2e-05  4e-05  0.9004  3e-02   3  3  3 |  0  0\n",
      "16  +1.916e+02  +1.916e+02  +4e-02  7e-07  8e-11  7e-06  1e-05  0.8174  1e-01   3  3  3 |  0  0\n",
      "17  +1.916e+02  +1.916e+02  +8e-03  1e-07  1e-11  1e-06  2e-06  0.8917  9e-02   3  2  2 |  0  0\n",
      "18  +1.916e+02  +1.916e+02  +2e-03  4e-08  4e-12  4e-07  5e-07  0.8588  2e-01   3  3  3 |  0  0\n",
      "19  +1.916e+02  +1.916e+02  +2e-04  3e-09  3e-13  3e-08  5e-08  0.9309  2e-02   3  2  2 |  0  0\n",
      "20  +1.916e+02  +1.916e+02  +2e-05  4e-10  4e-14  4e-09  6e-09  0.8768  1e-02   4  2  2 |  0  0\n",
      "21  +1.916e+02  +1.916e+02  +4e-06  6e-11  6e-15  6e-10  9e-10  0.9089  6e-02   4  2  2 |  0  0\n",
      "22  +1.916e+02  +1.916e+02  +1e-06  2e-11  2e-15  2e-10  2e-10  0.8362  1e-01   2  1  1 |  0  0\n",
      "\n",
      "OPTIMAL (within feastol=1.8e-11, reltol=5.1e-09, abstol=9.8e-07).\n",
      "Runtime: 6.538894 seconds.\n",
      "\n",
      "final objective value: 191.6347201927456\n",
      "CPU times: user 6.51 s, sys: 291 ms, total: 6.8 s\n",
      "Wall time: 7.5 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "import cvxpy as cp\n",
    "x = cp.Variable(shape=n)\n",
    "tau = 2*cp.log(1/p - 1)*sigma**2\n",
    "obj = cp.Minimize(cp.sum_squares(A*x - y) + tau*cp.sum(x))\n",
    "const = [0 <= x, x <= 1]\n",
    "cp.Problem(obj,const).solve(verbose=True)\n",
    "print(\"final objective value: {}\".format(obj.value))\n",
    "\n",
    "# relaxed ML estimate\n",
    "x_rml = np.array(x.value).flatten()\n",
    "\n",
    "# rounded solution\n",
    "x_rnd = (x_rml >= .5).astype(int)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Evaluation\n",
    "\n",
    "We define a function for computing the estimation errors, and a function for plotting $x$, the relaxed ML estimate, and the rounded solutions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib\n",
    "\n",
    "def errors(x_true, x, threshold=.5):\n",
    "    '''Return estimation errors.\n",
    "    \n",
    "    Return the true number of faults, the number of false positives, and the number of false negatives.\n",
    "    '''\n",
    "    n = len(x_true)\n",
    "    k = sum(x_true)\n",
    "    false_pos = sum(np.logical_and(x_true < threshold, x >= threshold))\n",
    "    false_neg = sum(np.logical_and(x_true >= threshold, x < threshold))\n",
    "    return (k, false_pos, false_neg)\n",
    "\n",
    "def plotXs(x_true, x_rml, x_rnd, filename=None):\n",
    "    '''Plot true, relaxed ML, and rounded solutions.'''\n",
    "    matplotlib.rcParams.update({'font.size': 14})\n",
    "    xs = [x_true, x_rml, x_rnd]\n",
    "    titles = ['x_true', 'x_rml', 'x_rnd']\n",
    "\n",
    "    n = len(x_true)\n",
    "    k = sum(x_true)\n",
    "\n",
    "    fig, ax = plt.subplots(1, 3, sharex=True, sharey=True, figsize=(12, 3))\n",
    "\n",
    "    for i,x in enumerate(xs):\n",
    "            ax[i].plot(range(n), x)\n",
    "            ax[i].set_title(titles[i])\n",
    "            ax[i].set_ylim([0,1])\n",
    "            \n",
    "    if filename:\n",
    "        fig.savefig(filename, bbox_inches='tight')\n",
    "        \n",
    "    return errors(x_true, x_rml,.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that out of 20 actual faults, the rounded solution gives perfect recovery with 0 false negatives and 0 false positives."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(20, 0, 0)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtEAAADaCAYAAABgkR7tAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJztnXm4HVWV9t9FQpgVkSA4QFpbFBEVie3QILRt0AZtUexG8RPT3TSN2No2jXZoUXACaZXRIAaQMCiKPEijIRLmhCQQE6aEMZCJhAw383hzp/X9UXXOPadOVZ1d866q9/c897n3Vu29a9WuvfZetWrVKlFVEEIIIYQQQszZpWgBCCGEEEIIKRs0ogkhhBBCCIkIjWhCCCGEEEIiQiOaEEIIIYSQiNCIJoQQQgghJCI0ogkhhBBCCIkIjWhCCCFtiMgFIsL8p4SUBBE5TkRURI4rWpY6QSO65ojIeSJyUtFyEEIIIYSUCRrR5DwANKIJIYQQQiJAI5oYIyJ7Fi0DIcQfcdi9aDkIIWZQZ8sPjegSIiJ7iMizIrJQRPZq2b6XiLzo7uuqmG7M424AvuTGUqmIPOjuG+/+/zcicpmIrAKwzd3nGy8ZFJMlImNF5I8islFEdojIHBH5RKJOIKSEpKW7bh0VkatF5B9E5CkAOwF8zrPvJBGZ36J373X3f1FEnhORXhGZLSJvz+J8CSk7Bejs34nI465uLhSRz/u08zYRuVtEtovIKhH5MYBRKZ0yicDIogUg0VHVHSJyGoBZAH4M4Cx3108AHALgQ6raa9DUFwH8EsAjACa521Z7ylwJYCOACwG8OqqsInIsgLsBzAfwAziTxikA7hSRz6rq7VHbJKSspKi7DY4B8FkAPwOwCsBzLfs+COBEAFcBGAAwAcBdIvI/7t+TAOwO4FwAvwHwnpinRUhlyVlnPwAnvPJqANcBOB3AzSLyhKo+CwAicgCAhwDsCeCnANYBOA3AuFgnSBJBI7qkqOqfReQiAN8WkTvczWcC+J6q/tmwjZtF5FoAi1T15oBi2wAcp6oDUWUUEQHwCziTz0dVdcjdPhHATDgTEo1oUivS0N0WDgPwXlV9ImDfYar6EgCIyGoAN8BZ/N+qquvc7X0ALhaRo1R1XoxTIqTS5KizhwN4V4vB/DsAywD8M4BvuGX+G8DrAHxYVWe45SYBeDKiHCQFGM5Rbr4P4DE4d6y/dP/+QcrHuCaOAe3ybgBvA3AzgP1EZH8R2R/AfgCmAniziBySkpyElIm0dHdWwGIMAA80DGiXR93fdzYMaM/2t8Q4PiF1IS+dfbbxj6quhuOpfnNLmU8AeKxhQLvltmP4aTLJEXqiS4yq9ovIeABPwXlce7yq9qd8mJe6FwnkUPf3dSFlDgCwNMExCCkdKepumH4u8/y/yf39csD218Q4PiG1ICed9VsLN8BxPDU4BP5PcJ+PIQtJCI3o8vMx9/dIOI+Jnkm5/R0+24I+wjDC83/jSccEAEGPian4pK6kobt++tlgMOJ2iXF8QupEUTrr1U2/NZj6WwA0okuMiLwDzuOk3wB4E4Cfi8gMVV0ToZk4XyXb4B5/X1Xd2LJ9jKdc4457i6reG+M4hFSSlHSXEJITFunsUgB+2XQO9dlGMoYx0SVFREYCuBHOm7lfAfAlOG/r/iJiU9sQ/THui+7vv/HIc6an3DwACwGcIyIdmT1EZHTE4xJSelLUXUJIDlims1MAvFdEjmmRb08AZxQgS+2hJ7q8fAvAUQBOVNX1ANaLyDcBTBSR01T1RsN25gL4qIicA2A5gDWqen+XOtMALAFwrZtfdgeAU+F5nKSqQyLyz3BS3D0jIr+Ecxd9EJz0W4fAeRuZkDqRlu4SQvLBJp29GMD/A/BHEbkCwynutuUoA3GhJ7qEiMiRcJT6WlW9q2XXzwHcA+ByEXmjYXP/Aeft/AsA3ALgO90quNk6TgLwrFvvm3AM5Qk+ZR8G8FcAHgbwbwAmwrljHgRwvqGMhFSClHWXEJIxtumsm7HjWABzAJwDJ8/7/XDWYZIzohonJJYQQgghhJD6Qk80IYQQQgghETGKiRaRD8N5bHAUgNcD+CdVndylzhFwPmv5VwDWwwnA/77S9Z05IjICQLeX9nao6qYuZQghOULdJaRcUGfrjemLhXsDWADn7dSuAfQi8io4sULTAbwPzlfrJsMJfP9pHEFJJN4EYHGXMjcAGJ+9KISQCFB3CSkX1NkaY2REu8H0dwGAiEw2qPIFOOlfvqSqOwAsEJHDAJwtIpfQG505qwCM61LmlTwEIYREgrpLSLmgztaYyC8WishWAP8eFs4hIjcCeK2qntiy7X1w3iZ9s6p2u2sjhBBCCCHEWrLKE30gnJzDraxu2ddmRIvIGXAThe+1115Hvf3tfh/jCWfdtj68snEH9ttrFN6w7x5t+5Zv2IEN2/vwxn33wGv2GuVbv7d/CAvXbMFuI3fBoa/bxxF4cy/WbNmJA/bZDWu27AQAHPGGjm+GYP6KTYH7wsotW78dm3b04+D99sSr99gVALBw9Vb0DrR/+fPAV+2O0fvsFto2AKzc1Iu1W3fiwFftjlWbe41kAoAVG3dg/bY+vGHfPbDfXqOMzycKYW1mcbw4PP3KZgyp4h0HvQojdhFs3NGPl9dvx6v32BUH77dn5PbSOq958+atVVWrPkyThs6SYTbv6MfS9dtxyGv3xKt237VocUhCbNPZNPQ16Rz+wuot2DkwhENftw92G7lLW72D99sTywLm2ua69urdMXrv8HXQK4efXI1trZjO0a3tRZnfTeRKwobtfVi+YQdes+covPE17fZPmG2UN2n2w/ptfVix0f+co5JEX7PyRE8D8LKq/kvLtkPgfKDjg6r6SFDdsWPH6ty5cyPJBAA3PbIU375jAb7w/oPxw08f0bbvnN89idvmLcf/fvZd+Mexb/Kt/8LqLTj+0ul46wF7456zjwUAXHbvC7js3oX42t++FVfctxAAsORHJ3bUHTNhSuC+sHJn/Woe7pq/ChNPfS9OfNdBAICPXzYdz63a0lbvf054O8748FtC2waAC+96FpOmL8K5f/d2XDT1OSOZAODc2+fjljnL8MNPvxNfeP8hxucThbA2szheHI44/25s2TmApy44Hq/afVfc+eQr+Notj+MT7zoIPzv1vZHbS+u8RGSeqo5N1EiGxNVZMsxPpz2PK+9/EWePOxRf+9u3Fi0OSYjNOhtXX5PO4R/5yYNYtHYb7vuvY/GW0Xu31bvy80fiq7c8jhPfdRAmeubaH/zxGVz78GJ864TD8K8ffnMkGf3kamxrxXSObm0vyvxuIlcSbp37Mr5521P47FFvxE/+4d1t+8Jso7xJsx9umbMM594+H59735vwo5PflUiuJPqaVYq7VXA8zq0c4P5eDUIIISQBa7fuxJgJU/DAc2uKFoUQUlOyMqJnAzhGRHZv2TYOTnD9koyOCQAI9asbON19i2T8HqR2ESzq4eNKy9c9Hbz9wG4hxD4aj4Inz1pSrCAkEqHrTMi+buskCe/bqvZe0XaLkREtInuLyHtE5D1unYPd/w92918kIve1VPk1gO0AJovIO0XkM3A+CZ1ZZg6JuS+sjBjVjI9f+yLxjxm3ZoJDVgsJ/ZeQzCl6QSAkM0Im1LA1iOtTd5LaP2XElvMy9USPBfC4+7MHgO+6f3/P3X8QgGbQrptUfBycD7PMBTARTn7oS1KRmhBCKoQtC4IJyzdsx7SnVxUtBiGEFI5pnugHETLPq+p4n23zAXw4rmCEEELs44TLZ2Bz70DhLwI3oPOeEFIUWcVEE0IIqSCbeweKFgFAubz3hJBqQiOaEEIIIYSQiNCIJoQQQgghJCKVM6LDU7x0j57zSx6Sdcxdtzfyox4/7hv+jC108XYEO4YQQlIieEINW6OZuaY7dey/olMfVsaITpoix69M5ql1/I6ZcnsZVqsc3n5gaiWSN0UvCIRkRXgatuC9SdK+1oWwPqpq99lyXpUxogkhpLTYsiIQQggxhkY0IaQQxl8/BxMfeDF2/WP+934c9+MHUpSIlJGMvt9FCCFdoRFNCCmEB5/vwY/vfj52/ZfX78CSddtTlIiUCT7mJ6QalPlGmEY0IYSUjDETpuB/fj+/aDEqxR2Pr8DJP59VtBiEkBJBI5oQQkrIrx9dVrQIleLrv30C85ZuKFoMQkiJqKARnSzFi1+RrJ80lPdBRjVhhgRCCLELzspJsbcHk9hYRUeCVMaIDk2RY5TErbNMARnufF/SNx0kjfOMagQ2j1n0aCyYoBhLGtX15OX12/GpiTOxcXtfbsesuQqSGuA3xhtTr+++bMWpBM0+8u2/avagLedVGSOakLSxRUlJMUx84EU8+fJGTF2wqrntivsW4tmVm1M/FkcaqTp8EZRUERrRhJBacdav5mHMhCmR6w0OKS655wV8auLMDKQihBBSNmhElwCGExQLH7FXi7vmr+peKISBwaGUJCFJoF+TkGpQ5iWWRjQhhITAmyhCCCF+0IguAYzNLRaG8tWToOtOm5oQQghQQSM6zGtktPj5FMo6nKLMX+upIrwc5WX+8k1Ys6U3k7Z5L0VIcXBeTobN/ZfEBir6tCpjRId5C008iX5lsvZA+r2t7JviznCYhKUJCq3XTI1XbwI9j3XvmBLxyZ89jI/+9KGixYgNhxqpOn5jPDTLKu9eu9Jc+0P2VQ5LzqsyRjQhaVPZyafibO4dyKTdLA1cjjVSdTjESRUZWbQAhBBSdv60YCWWrNuOA/bZDe98w6tx6Ov2KVqk2sAnRYSQojA2okXkLADfAHAQgKcBfF1VZ4SUPxXANwEcCmAzgHsBnKOqyfJL1RAuEsXC/q83JuFUZ978WNv/S350YlbiEBd67wmpBmVeYo3COUTkFACXA7gQwJEAZgGYKiIHB5T/awA3AbgBwOEATgLwDgC/SkFmQgjJgYDPwJd5xieEEJIapjHRZwOYrKrXqOqzqvpVACsBfDmg/AcBLFfVS1V1sao+AuBKAO9PLnL9oMelWNj/hBBCCPHS1YgWkVEAjgIwzbNrGoAPBVSbCeAgEfmkOOwP4HMA7koirAmhKe4MPEh+Reh5qhe83IQQYhf8cm8ybLZjkshW9HmZeKL3BzACwGrP9tUADvSroKqzAXweTvhGH4AeOM9Gv+RXXkTOEJG5IjK3p6fHUHRPG2H7TFLc+csVSxZTfI/ps9V0kMSVNm5qvKoR1H917xc/0tDZspFLPncONpIBNumr3xAfTtHWuZMfG+vO8Bru13/VxJbzipLiznt1xGebs0PkHQCuAPB9OF7sj8MxuH/h27DqJFUdq6pjR48eHUEkQrLDFiW1EepsutBQSIfe/kH0Dw4VLYZ12KCvDIsjVcQkO8daAIPo9DofgE7vdINzAcxR1R+7/z8lItsAzBCRb6nqy7GkJYSQnKFz2G5avZdv//af8JcH7I17zz62QIkIIXWhqydaVfsAzAMwzrNrHJwsHX7sCcfwbqXxP+9HCSHWQ8+Z3QR5719cszVnSYgJvBklQZQ53t00nOMSAONF5HQROUxELgfwegBXA4CI3CgiN7aU/wOAT4nIl0XkzW7KuysAPKaqy9I8AZIducSAloAyK3hdUdVcHusPDim+/pvHMz9ON1Zt6s39mJwfCCF1x8iIVtXfAvg6gPMAPAHgaAAnqOpSt8jB7k+j/GQ4afH+HcACALcBWAjgU2kJTgghQZx/59N467emZn6cxWu34o4nXsn8OGHMWbweH7joPtzx+IrU2pz54lqce/v81NojhE92SBUxfrFQVa9S1TGqupuqHqWq01v2Haeqx3nKX6mqh6vqnqp6kKqeqqrLU5TdX84Qr6GJR9HPu5K1v4UOHbvwjgG+9FU+bpy9tHuhiNiqps+t2gwAmLd0Q2ptfuHaR3HLnHo+NKSH3U54WZJh8xPVRCnuCj6vKNk5rCb8Lre7EZR1Ojv/Y5pti9pe1EWgcUhvvbpNWkFjoGglJfWBI41UHf/51Jl7w9LfkWAajh7fnq1o/xVhs/lRGSOakLSxREdJArL3KqYzSMo41up2k02SwSd6pIrQiCaEEB+aT2gKlYJ0I21jnjcHhBBTaEQTQggpHWX03tcZhsWRKlILI7rsk21RL7pwynOgZ6q88NplB7u2nvDFS0KGqYURTZ0nhCSF8wgh8WFMNKkilTOiwxY6k0XQr0jWi2e3x1y2vIVaF7xXg91fXmj3kqhwzJAqYrMTIJFsBZ9X5YzoKhL18VncAemtxsd2DuwGEkTaN1gca8R2ko7RUEdXsqZrD+eP/KmMER32qMhkofMrkrUH0veYidqLV5uebgf2AvElh5WpjGPPlptsS8QgCQhbgsqoG3kT3n/V7EFbzqoyRjQh6WOLmpK4JDH0eG9pN1ldHltuDqoG9YlUERrRhBBCCMmUzN8t4s1PaSlz+sORRQuQB2W/AS5qbijvsE4X9gPJir/64b344gcOKVqMWFAv6gmvOyHD1MITTaUnpJ6kqft+3pKkN+hrtuzET+95IWErJE24XmQDwzlIFamcER02AZpMjn5e36wfNXTzNHPyyRfv9WD3l4+8dIYGFyH5wHCNZNjce0kubdHnVTkjmhBCGqSx7uY5SZcpNpA2DSGk7lTHiE6YIsfPc5V1ahjf1HI+20wXq0bVuGub9zh1WySDvJc16wbiYqL/ac0Qtj5tKoP3L+0bjxKccqEkHRN+1SVsn6W6YSO+V6ai/WfLuKiOEU1IytiipIQQH6ifhJCCoRFNiIcyPVIn4fBaJoNeWVIWOFbLS5kvXS2M6LJ7FIsaYDRAHDg558+YCVNwwZ1Px66fhcr7P2ou+eSSAM4P9YRXnZBhamFE0wgipHxMnrWkaBFIDeHNASHElMoZ0aEGs4E17TuBZv6lpWzbJxHh9agM1K1ksPuILVCXk2Fz/yV5WbXol58rZ0RXkfo+MLaDGj+xrzXNbDc5ztE2L3ReyiQrIYRkgbERLSJnichiEekVkXkickyX8qNE5HtunZ0iskxEvpZc5IDjhcpiUr+zUNbGk1/zfttM16qwNEGh9QJS49VvkfS/4PXrB5I3UWKri/a82Eba3cHuDSdp//h++bM5/v2+CkovRjca/ec3N9jce0lsLFucWyNNConIKQAuB3AWgIfd31NF5B2quiyg2i0A3gTgDAALAbwOwB6JJSYkJyzRUWIxVR8jzqJs51nSuCoXdX4Jl1QXIyMawNkAJqvqNe7/XxWRjwP4MoBzvYVF5HgAHwXwFlVd625eklBWQnKBjqjyIyK5uBTTPgLHHqkqWT9Boe6UlzJfu67hHCIyCsBRAKZ5dk0D8KGAaicB+DOAs0VkuYgsFJErRGTvRNLGpOweCz5eLBpegLJC3UkGu494YfYSQoYxiYneH8AIAKs921cDODCgzpsBHA3g3QBOBvDvAD4OYLJfYRE5Q0Tmisjcnp4eA5EIIUViu87mdduc9nHKfbtPbMV2fSWkrETJzuG9/RSfba3tKoBTVfVRVb0bjiF9soi8rqNh1UmqOlZVx44ePTqCSH5CdorU2GZy/+znucr6vpv39XZB72V30tTZLEjjEg6/qGvXgLBMHFICbNDXpDHRHPbJsPkJQpI5reizMjGi1wIYRKfX+QB0eqcbrASwQlU3tWx71v19cCQJiTVvodYXXoCyYvPC4Ydt0oYtblU15qt6XoSQ9OlqRKtqH4B5AMZ5do0DMCug2kwAr/fEQB/q/l4aVUgTwu5yTWKi/apnbjr5HdNnm/Gk3khzE3EpbvSP1+NWt8UkeAjVrCMqAG976gO1M18Sp7jzqR+WnpVOpO4EJwi0OytKEslsedfNNJzjEgDjReR0ETlMRC4H8HoAVwOAiNwoIje2lP81gHUArheRw0Xkr+GkyLtNVdekKD8hmWHz5EPyp243lbZD9SwXvFykihiluFPV34rIawGcB+AgAAsAnKCqDa/ywZ7yW0XkowCuhJOlYwOAOwBMSEtwQrKibCEAJJisDd+0DTnbDPUwXaCekChk/m6RbcpDzCnxpTPNEw1VvQrAVQH7jvPZ9jyA42NLliJ181jEfcxBzyshpO7w5iBbwpYZLkHJYPflT5TsHKQgok7qcReBjphoLiYA7PMOEnPSuHR5Xn4TI4LDkZSZOr6smhfsvvypnhHtl6JO23+HVvetn/GXljhzWAUvBwH4ZAagLpD0SKpNHIsJsbj/kjjsih4X1TOiCSG1p6z2b9ELQhTKJGsUqnpehJD0qYwRHbZmxl1Qs16I/WKXk6V8cYi6CASdZ90Wk5LaXSSENJ/yZBneVFaj3wpqNk8VTRbrQmP8+6ZoS/9wlUNCctzZ3H9JnvbZMmdWxogmJG0s0VESg7xyiNqSq7SOsOfLhS1GDyFpQiOaEA90bJWfsr4UWya5yyNpNKp6XkWT9ZNNXrfyUuan3rUwout2A8w7fkIcSjw3B5Lni8g2L24Wi0ZiwqUrIezA3KmFEV12oi5kreWTLLhcpBzYD+UjzTCLPA1JhocQ28nyaQnn2oSwA3OHRjQhAdCjT4i9ZKWeTDmaDZxPSRWpnBEdNv3FnRwZy1UvGuOEa2n5Kd1nvy2bDUI/+00FITlim26UDZt7L4lsRZ9XZYzoLO5ys3606idzGilfWgeVyToXlBqvboskPSUVIoNrmaU6MIwjPjSu8iWpHoR/sdAvRxt1oxuN+cNPF2zuPotFM6YyRjQhhHRA+yoRNt9H84uShJCioRHtYvNiQfKlMRa4RpPcMZiHbJmqbJEjbap6XlWHa3h5KfNT71oY0XXzWMQ925p1E7EUWybUsK+okeKxZZyQaIStM3Vbq9OGYWH5UwsjuuxEXSo04O+u9bwx0RGPW1W4WJeX0sXLWrYGlqz3SA4kHRN+02ljG+faZJRuvqsANKIJCYBOkfLCS5c9Rds7WXktiz6vqsL5lFSRyhnRYXeyJnOjX/3M51RO2lbCxbS8ZHHpchkOHHOE+ML5OBk2e/kTpbgr+LwqY0RnkuIu4ztnv+ZNt/m3F09geggcgvqPcXr5kuacaPG60UaUIZbnORW9QJlQAhFJF0LjpPMTo7SEx5nnJ0dUkshmy7pcGSO6ykReI1pWlSiLoDeeiouTQxkMCdKOHdMryZLM1lCqeyhJ50O/uN1mTHSilgmXqvypjBFtooShSd7dmn4fKslqYGrzt3Zsi9eenV9kLAvNMcAUdwTFvOlumyqGymObsKS0hK7bBmWc/fUdkOEfsHF/5yNKJJLYWM0vC6coTxwqY0QTQoiXNCbYMM8bb7SKgzf/hJCiqYwR3VjMwtY0kwWvtUizzYwWSmn+lo5tbRiuFn6eM5OagedXs0UqyPNYs26oBGU1bksqdqWos0fThLx7Rzy/SScmMdE29l8aMdFFn5exES0iZ4nIYhHpFZF5InKMYb2jRWRARBbEF5MQUhfSXKTTiGdvTNa+mXsqbm+ZhMAVRVlvlAgh1cHIiBaRUwBcDuBCAEcCmAVgqogc3KXeawDcCOC+hHIaEzcmOqx+1gtl0YsRaYfXg7SS54ulRmk4OT5JCUn6jgFHfTJs7r8kc1rR52XqiT4bwGRVvUZVn1XVrwJYCeDLXepdB+AGALMTyGhEGT936edJSfZ4I2698vVdFrAbSFFw6MWn6EWUJCdpGGbdCe0/i2eXJLLZclZdjWgRGQXgKADTPLumAfhQSL2zABwI4AdJBCTR8cswYlSv47PfXJ5IueEITojFHZhZhjuLz9kGkvaP72e/U2q77rD/8sfEE70/gBEAVnu2r4ZjJHcgIkcAOB/AF1R1sNsBROQMEZkrInN7enoMROrEm57Mv0zIvpBUK1kZk37HTKIEcesyD7KD93rYfAdfNGnobBCpxDGncO2K8ICVSRdLJGrtyVJf0yD0S8OGFnadx2N4isDutlFRNG2rOCnuOv4ohijZObyiis82iMhuAH4D4BxVXWzUsOokVR2rqmNHjx4dQSRCSBFQZ0nRWGgTWAv1lZBsGGlQZi2AQXR6nQ9Ap3caAA4C8A4A14vI9e62XQCIiAwAOEFVvaEhiWl4nsLijE0cS771M/JI+qXQ8zu+6WLhX7d77aCYaBvvXLMkyPNYt36oEql+Qjy9pgKx7f2EOoZ01e+MI5JzBzVVwjLdsAmTmGgbu0/8fbGGdb1/FENXT7Sq9gGYB2CcZ9c4OFk6vKwAcASA97T8XA3gRfdvvzqEEEIMsOWmrmgxLLQJCCE1w8QTDQCXALhJROYAmAngTACvh2McQ0RuBABVPU1V+wG05YQWkTUAdqpq5rmi48ZEh9XP2htjy6JIHJqXg6t06SmbJ9W2mGjLxCElJrEnlIMxETZ3XyLRCj4vIyNaVX8rIq8FcB6ccI0FcMIylrpFQvNF50EWjyqyfrTqFyZSRMoX2ooOHf1g8aRTZdLodhPVVVVrwicsEaOUpH3jYduNTB0I00O+4N0dky8WWkkBKX3TxtQTDVW9CsBVAfuO61L3AgAXRJCLJCButg9vUS4lpM405mjaVP4UbWzasojWjaRPd/w/aKaB+2yk6LEfRNmevFWBKNk5rKaZniy0TEganebv4TJNxc5oXDZTz/hsi9devvWqxnAmpWaOO1J2sv7aaMrt27Y2WyYOqShGKdqyPL5tiheR0DDW+FnksicktXDXqs3zKvbMKmNEE0JIFKJO3CVfZyNRBi9vds4NQggxozJGdFgau+EUL91Xhtb4q0b5rBaUplwBx4/envs7YhPi+V1X6n7+tlAnY9VLGYxXQpLgv0aH7etcJ0k7w/aPz3tWEVL85k4KMdFFx8xXxoiuMnl5zDpioutszbTAXigvvHbJMAmBKwrecBRD4s9+h2wrekyZYuvSWHRoQx2pnBEdNyZ6uH5nmawVhsPeLhrXg2t0eTG5djbqna2LMyFJSZzhLhUpAtqugd7ZfIrJ3gVjTLS1ZG1E+XpSikj5QmsRQGe4j82TDskeozR5KY2Soh9JBlEGHUhbxjoYVLYR+sU9O1XDMsrZSUmktmXOpBFdQdozjESo5ynLtYSUnTQMoqI9HbZSvLFpxyJaN5Jedr8nwh2ZkSzHVilL0n2VonJGdFwl9Etnl/V4bKZoac9xl7gN0/VVAAAcQklEQVS96BXjH7NKZJ3SkJiRl9Fq44Id96uqaUPzlOSJUYq2LuM+iV7YNxNExSBU1eL5LlaKO9ixXlfOiCYkLWhIlJ80DPLmY0P71qDMKXqBCicb4fjkgRBiSuWMaL80dlFiqlrL5hUT3SZfESlf/OSoIbZ8ApokJ81rmadRxRFYHFT/fAhbo8P3ZSlV2emejtfG9S2NmOiiT6tyRnQVibqIt8VER6jbERNNhwwAOx+DkeTwqiajeI+tfUZBHUg6H/rGRGvwPhuxVU5Lxao0lTOikyt4SoJYfkwShnNBbLxzJ9FIolsm1z/1z36n21xiijeUC6CGp5wLCefTLNdJW43iNLH5HJNIVvRpVcaIzsLgydqG8ms/0eONmALbkirGNmyedKqMbd2ey4t8lqlgmW4gbRsvJDqhYQj5iVFayhbG0SCJaLacVmWMaDJMW4aRKCnuPPeDtfRCEWI5uepl6NercpPCF1sW0bqROMWd37aSLTW2imurXFWmMkZ0Mz1ZaJmQfZ7freWz8kg222+LYU7SXsz0flS9Nso2oZNg4s4H/m1lPzA49oqHl6AIgnt9eG0OvzLJvnpXbsJTBHa3jYoiiY1lmvowaypjRBOSNjY/BiPh8MpVn6IXT0IIqYwR3TB4wr6kbWITtWWbC0m7kwZ+KenSiImOKm4zVQxND1Ixso5rT7t12+7bLI7mSB3Lur6y+PdzcO83129eoUCGbQm/fcG2UdGkERNd9JxZGSO6ykR+7Bzzq4sdj8OqtkoS0gLDmMpN0YtnXUl6X+pfvVy6aOtTEEvFqjSVM6KTxkD6v/SQ7tD0tsfF3C54NQgw7LnJajyE5csNr5eBMB7qbJ/aaiDVnSzXyTpcc5tPMdEn2xkTnQ5hk35cj0X2C0myryumJkWdV8wWvP3Abik/RU+wZacM/VcCEQm6zach4RyciLsSav/kJkV0koTo2HJelTGiwyjDQkDsg8OmxBjMsFHnhS29A/FkqShFz6u2LKKEkPpibESLyFkislhEekVknogcE1L2MyIyTUR6RGSLiDwqIn+fjsj+aMcffmVC0uj4pEvxS3vXWS9+apZGy3HzOge233as7g0GpYqpmxFpS8qcupNF/z+/agtuf2x5rLoNY+2UX8zu2Jd6qJdlWmebPHlQx3OOQtL+8VeZzrXQW76bqiUKCYh4TraNkaat4td/HX/YQxppCYu+FkZGtIicAuByABcCOBLALABTReTggCrHArgfwIlu+bsA/D7M8CbENujpqg4fu2w6zr71yVh1G1P0K5t60xOIJCb9zCjU+KKg44KUlZGG5c4GMFlVr3H//6qIfBzAlwGc6y2sqv/h2fRdETkRwEkAZsQVNgzp+KNln08qucB2WoqI53f6dKakS2Mej9qELaliiqbu50/yxc9wsCWNl4kuFO0BIuXEb2yFrUFcn7oTlsY3zDYqmjRiooueM7t6okVkFICjAEzz7JoG4EMRjrUPgA0RyhOXqEuVBvzdtZ43nINrJAD2g+1843dPord/MP2GLVx0gHyfyto89rO6PDafsxUkTXGXIGTDFmyV01a5qoxJOMf+AEYAWO3ZvhrAgSYHEZGvAHgjgJsC9p8hInNFZG5PT49Jk8HEjIlulgmLKUoJb3sc+HZBD1t3UtVZD1H7/3fzlmPqgpX+baXxjkGOw4Fjj2RBlvqaF1nqYV46nvXHn8IPXtyhu5FnPHvaRMnO4ZVUfLZ1ICInA/gxgC+o6lLfhlUnqepYVR07evToCCK1HidUhphtZuuK8n/0kuDxhrWp/MpBR4o7dkwgaehslphcurRf4iUFQS9EV2zQ1/A1OmwfJ+JuhPWRzd2XxhcLi8bEiF4LYBCdXucD0OmdbsM1oG8CcJqq3hlLwhQo9O6PlBYOm/IQdK3S8FLkMVnbaigk/XhVGanoaRFCMqCrEa2qfQDmARjn2TUOTpYOX0TkHwHcDGC8qt6WREgTmjFVBmnsAvZ6fg8b32H14iwk3nbbY5gTpHzxS9MX8yuNSWUpI2WLyyPDeK9ZmkZpVuPBr1lbxl7YEzFL7X0jbn9sOa6dsahoMUpN8iHa2ULY+t1cLyO3ah+ZzSXaab94j2njep5ENlvWa9PsHJcAuElE5gCYCeBMAK8HcDUAiMiNAKCqp7n/fw6OB/ocANNFpOHF7lPV9emJT0h2lNlYqBsKYPoLw7GeRvnRLVxUiDlxrl4jzeHpx7y5Yx/VvTioiaSsGBnRqvpbEXktgPMAHARgAYATWmKcvfmiz3Tbvsz9afAQgOOSCBxEWBq7hlfKzCgaLhStXnSGZW49evKY6KjycvFwoNFsB3E9C5Omd3oYk3gpikidZJRaLkfXS9ixaPiQePis0c3fIet3liKVnOEnb37919hjXw+mERNd9Lpt/GKhql6lqmNUdTdVPUpVp7fsO05Vj/P8Lz4/x/m1TcKJuma2lU/y1itXSQD0WJYBVfXkW482sy7q2YoxE6bg6Vc2dbbN618r+A5NOMm7xyfkoPm7HH1v7xCxVrDKEiU7RymIHxM93EK8euZ0pLhLt3mSEHsnSBJEnEvWep3vecZ5R/qOx1fEbqOVm2YvMaib/VxjyvIN23HDrCUd26kLJApZ3oBkORRzS3GXz2ECjm2vMidKcVfwaVXOiCYkLWx8/EWC8fM+m86vjappfUXw2//3dOQ6abNgxSb0DQwZlT3tujk4/86nsX5bn3H7tnhs7ZCCdIOzKakilTGiszB4so618Ws+UYxQzD4oOqbIFrz9Z/Ode5WJ1esd2TkCigUYfo1r75sxI4dxEFcFr5+5GL94qDMWfNm67fjElQ/je380M+Y37egHAAwOaVeBqj5dWHJvUCvKNqZsWxvC+s9mZ1Ayyew4r8oY0Y1BHZqOLqx+SHq40BARQ/l8j+XdgISPNRp90JqmL8JXGr0Ghl3TRPb49R8pBwptm1ID80a36nfL9jBPdFo8vmwDxkyYgjmLkyUoahXxu394Bis27ugos2G741F+anlnjLcfzfOvyNif+eLayHXoTDAjbsikhpQJ22fStrM/QXrYiOM+7qGyenoT3rfdbaOiCE7MZ1676NOqjBFNCKk3fkZQ1BvD1sU0baOqYdg9+PyadBtOhaYVnSpzFq/HmAlTsGZLb7oNI9wo+MK1j6Z+PJIdNhp4hJhQGSO68cjC/1Pa7b9D25HOv7N6HNJsteWgaYRzRJV3OFVMvV0x3n6z+TEYaUfVkyoy4NINBYVzeCqoKq57eHGz7ZAjG8vYOEYy70v6LN+wPda8Y2L4/NLtw3lLNkQ/AKkUoSEHBunKar48hdK0cXztn2xT9SYhjXCOok+rMkY0IaS+KMxeLGwP1wrOjrFgxWbjl/KiYpPX7Z5nVuPoix9Az5adAML7y8sDBXnU12/rw5gJU/D7x5cXcnxCCGlQOSM6bkx0WP204wQ72uuyqkaO12r7dHmEet5H3zat9jniPe2adkNhxBl3HZ7ooHIBuuQt3z+UvgFtY9zxghXtMdONrg/z7jSk/05BGUiWrNsGALj5kWWFHL/uJF5jQ9pMEhOdhMjfYoh7nJj1jNsvqP+SkiyevVgqZ0QTkhY2Pv4iZgRNrEFz9fCLhU4Bk9CQqEhGccdJ8J6b18AfUsX46+ckjuP+4/yVOOr796B/MPnNifdypJ7H36LrUyXqHi5IqklljOhQ/Yypu1mrvG/8UgqfwYxcr/CoIjvgHF9eFBrwYmHw/2G2UuuCn5ZRJSE2dJSxFyTPph39eN8P78Xjy8zjjzvSOnra7h8cwoPP9+Bfb5wLAFi7dWes/pjy1Eqs29bXTKWXBBpj1YOx0Mkoa/8l0WVbzqsyRnTo4yCTR0We323bQh+RxHv03PbbZ5+vcIbtwtBQGC7jnwKnbh6ZoPOvWz+UEecatbygG1QuQCN28bz0t4ungc29/sZflLHRaLIxZ5ik7No5MIhFPVvb47MD6s1bugE9W3biyvtfNJepwxPtv78h6/l3JgvhCHqx04SdA4P42KXT8ciidYlkCILOBDPCrmDYehg27sPmWpNUs93k6kbUunHDD7JaS/xsiY59Fq5jzTERq27776KojBFNCKk3Jp6JbuEcDSPPa1Ct2bwziWhtxzCd9H/+0Et423l/wkd++lCb8drNmIhiqHpvFjrbam9zcDDhipWg+rJ12/H86i340dTnksnQwprNvcae+97+Qdz37GrfvNwkGTa9J2BCuaQlWVIZIzo0RY5B+hxP0ba/s3psMJxCr3NboFAG7UV1qNicAidPvOdvknaJmLO9bwB3zV/ZtVwsrwT8h/1HL3mo7f9L73lhuE7LgcSzzXvNuxmbJoR9FdGP2x9b0fx79kvDHw8ZCmhgF3c2j+KZ+b8nXmn73+thG/J4ipLqQhLjIws9PP6y6fj0VbPatgUZdBu29+FfbpiLGS/0pC9IDQh7dB+aojbjVLNVwM+W6NhnYfelEc5R9HlVxogmhNjNBXc+jbN+9RieeHmj7/4f/PEZfPG6R3H3glXRG/cYf62Tc2sqtmvdvMXd6LyhSj5TN5q4zlCGVsO9PeTL38h7eb3jIX3ohR5MfMAspGPhmq1t/3tvIryPTHdJ2A9JwjmyiIXeuD04RntRz1b09g82/x8KuMEihNSXyhnRSVO8hMUUZUXX5iMHbLX8GUH4omOLbKHuseFZ8cpG56t1WwLii699eDFmLFyLb9z2VHPbrXNfxuHf+RPGTJiCeUuDP5f9oscYbOWXJkarhHuJvZ7ohau3dG8zBq362up5azU+gzzRrUx7ZnWs43uNXO//iT3RUVOJqWLa06swOKS5+iF39A3iIz99CP/52yfaZAH4YmNo3HPM+iafps4y5CNqjHPsz35nHAgSHq+e6aETkUS2os+rckY0SU7N14gm7IZ0aYyrQR8rcMpT/mEe37ztKWzrc7yBd3pCD1q5YfbSwHEbaPSEhHN4Pa7eR8njLp0eKEsQHV9FDFjyVBULVmzC8y2G+lDbe4XZrRreS+NdoJIakFElnzJ/Jc64aR6ue3hRoBc8bn88t2qzf3s6nCf84YVr27YDyb3xdYW9RqpIZYzo8Ax3cXO/Zav2fnIlifuKW5NrggO7IVuan7322DyDQ4qv/PqxrvVNPLB+dItn7hsY6vjoSFB8vJcoIplm4bz50WX4xJUPB5aN43nZtL0f46+f0+ax/9n9CzvKdfNEJ40NH4p4Ede6X1JcvmFH4DVYsGIzxkyYElmW3v7hO5OZL65FX0sO68ahtuwcwOrNzhOURl+kER9P2glb9xgL3Z3QFHf5iRGZJLLZcl6VMaLDKNubv14SRHNEqustW/RjkqLR5u+ad0RKjGjmSW7vTz/PtB/d4mmDFttunsMfTnkGv/nzy2hIF9aWl1kvrg3dv2l7P7btHHDabGnyK79+DEd+7x7fOhfd9WzHtiFV3PrnlzFmwpTmJ7qj8MuZi/Hg8z1tL1r+ZNoLHeVU2x9tey9N0oUr6pwywrVYh1S7Xsd1W6P1S2trv360/euHrWK+/8L7XBnceras3gWRdF3wDZlshHOUZa61VMy6fmW4SGphRJsQP8Ypv2ORfCjNRF4yGp5o7xe1TV8262Zrtxo3fll2vCgUyzdsx4Mt2Rb8snOELUzXz1oSuE9V8e7vTcPRF9/fIceUp1Zie9/wS2s7Wryirdtbt/12rmPoL1q7LfCYQYzwuE+vn+kfJ37vs6vxF+fe1fTSes/da8h+87YnMWn6S8ZyDAb05QmXz8A10xd1bG+OGYMh8jc/edBYDsBzLp6XONXnw4ra9ETX3IoOIfY6msKUmyiuNnL5eAfj2h9EeTumFkZ02R8HRZU+dlhHzHpVZThjIHsmDRp2XLeQgSC27RzA5ybNxks9/i8RRo6JBnD0xQ9g6brtbbK8sHpL+8dN4L/47RzoNHZb6XdzKm9wM0CEyXHxn8JzH2/a0Y8Rbv1B712IAa029NadA/juH57xLTfVkxmlwxPtOYdb5y7HhXeZ520OeurwzMrN+KGPB775ERzVruNkc++AsRxT56/EJ3/WEjLjabrfp4+HPdGcD7KCc20yODbzpxZGNCGkeBoT/PptffjW7+c3t4dl1mhlxsIePLJoPS4yMNrWbeuLJWPPlp04/tLp+NYdC5rbRMTXgDv71idDPUuX3tseLpF0fWvkgd62M9x4969rdnBvqVZPdP/gEJasi+4Fb6V/cAg7PJ72sDCMhtiDQ5qqF++CP7R/edHrWbxmRqdXvJmdIz0xCCElZ2TRAqRN2GMWk0cw/ul30sUrR9fP/0Zu37ztsOPUNayhI8VdTfshbRoG0YTb57dt//ufzTSq3whJCPLEBnuxgrJgdG5rGKhPtuSyVvUfAVOeWok3vmaPQHl//mB7mIOJ8RU20hpe2fNaDPxubHW9s60hCGEx6GGfAf/hlGcxb6nZ1/2COP2GuVixcQeW/OjE5rYr7ht+wbG3fxC9/YPYd89RbXIPafQc04ND2hbG8mjL58JH7tLuP2ptWlWxxCdkplGE4RzBmK2xfvUM6mc4DUdPvZiNHEkJTz9oqdBI2p/FnpexJ1pEzhKRxSLSKyLzROSYLuWPdcv1isgiETkzubiEkLKS/FGtU38gYoaHTTuCP6jRcQQfEX/+4Ev49FX+hr7pS5GBjUfAG9dsQiN+ekTLscMyZDy+rP1DOK3n9+ji4DzdXp54eSPuf64zV7XfJ7NHjhhehk75xWy853v34EdTn8PA4FDTgz4UcCMTxkkTZ7aF3Jwy6ZHm3x4bumMR9zOUmZ2DEOLFyIgWkVMAXA7gQgBHApgFYKqIHBxQ/i8A3OWWOxLARQCuFJGT0xDa/5jx9oW2Ga+aefs+B0iyzqZ1nnVdIxhPli1ewyVyfffyDAwGmFMBl+/PS8y9p77e6b7BwC/brdzUa9TutKdXJdarwPM2oDXe2C/eN4hrZwy/gLjrCPMzOGniTPzz5LmB+2e9tBZrtjh9N7LFKn1yuZNq8OqHXsKDz/c0r7mTNcT48ACA+Ss2tcnfSocn2mOi+4W/NLqN00Q8Yq9P7O+uhKYItLn/Etk7dpyY6bJ2NoDJqnqNqj6rql8FsBLAlwPKnwngFVX9qlv+GgA3ADgnuciEkDKSdNJreGJnL1qHH/yx88W4NKbUOUvMva1RmDJ/ZeLFbHZLOEISfna/2SfBAeDhF4czl+w6IvpdUNA5n3rNozjJDeMJitfe3j/YluIuTvquW+Ys893ezZvs54lueNFtWbwJIcUj3SYmERkFYDuAz6vq71q2TwTwTlU91qfOdADzVfUrLdv+AcCvAeypqoHPV8eOHatz5wZ7MK6ZvgiX3duZ33RH/2Dz7em9Ro1o27et5UUW774Gg6rNtE6NMtt8Uk156yuGU1IFte2VY89RIyABcvkd06TtoLp77Dqi64LhlaPxv0ldU7aF9FHYvjxpyLH7rrtghAgGhhQ7B9rHRJz2utWdde7f4tV77Bq4X0TmqerYyALkRJjOLlm7DSdeMQNA8Nguij1HjfBNJ1c1WnU6y2O0Ynq8KLLtvusubR9IiSqb6XH22HUEdvQHl530xaNw/OEHhrZhs852W2M/c9VMPL+q87P2YeuCd23zo1Fmt5G7NJ8+NLaNGrlLMyNOnPUb8F+L/eTyGwdhcvvJ0jp3RFmbvWt/Wmte/6A2PxbkbTPMNsob73l719wo7BwYaob2hZ3Xh/5yf1xzWrgqJtFXkxcL9wcwAoA3wG01gI8G1DkQwL0+5Ue67bV941dEzgBwhvvvVhF53kCm8K8ckDRgP2fP/vt+v2sfH5KLJBGgzloJ+zgHPnaxUT9bpbPUVythH+fAM8D+134pO32Nkp3D67IWn23dyvtth6pOAjDJVBARmWvrXX6VYD9nT1n7mDprH+zjfChjP1Nf7YN9nA9Z97NJkNtaAINwvMutHIBO73SDVQHlBwCkE9hHCCGEEEJIQXQ1olW1D8A8AOM8u8bByb7hx2x0hnqMAzA3LB6aEEIIIYSQMmD6uvUlAMaLyOkicpiIXA7g9QCuBgARuVFEbmwpfzWAN4rIZW750wGMB/CTlOQ2fixFEsF+zp669HFdzrNI2Mf5UId+rsM5Fg37OB8y7eeu2TmaBUXOAvBNAAcBWADgP1V1urvvQQBQ1eNayh8L4FIAhwN4BcDFqnp1irITQgghhBBSCMZGNCGEEEIIIcQh4TfECCGEEEIIqR+lM6JF5CwRWSwivSIyT0SOKVomWxGRD4vInSKyQkRURMZ79ouIXCAir4jIDhF5UEQO95R5jYjcJCKb3J+bRGRfT5kjROQht40VIvIdqclnvUTkXBH5s4hsFpEeEfmDiLzTU6a2/Ux9NYf6mj3U13Cor+ZQX/PBep1V93OqZfgBcAqAfgD/CuAwAFcC2Arg4KJls/EHwAkALgTwWThfnRzv2f/fALYAOBnAOwHcCid+fZ+WMlMBPA3gQwA+6P79h5b9r4KT0vBWt42T3Tb/q+jzz6mP7wbwT+65HwHg925/7Ff3fqa+Ru4v6mv2fUx9De4b6mu0/qK+5tPPVuts4R0UsTMfBXCNZ9tCABcVLZvtP+5kOL7lf4Hz5chvtWzbwx00/+b+fxicj+P8dUuZo91tb3P//zKAzQD2aClzHoAVcGPu6/QDYG84edU/Wfd+pr4m6jvqaz79TH0dlo/6Gr/vqK/59bVVOluacA4RGQXgKADTPLumwbmzINH4CzgfxGn2p6ruADAdw/35QTiTQ2s+8JkAtnnKzHDrNrgbTgrEMVkIbjn7wAmT2uD+X8t+pr6mTi3HUQ5QX0F9zYBajqOcsEpnS2NEw/nO/Ah0fiVxNTq/jki60+izsP48EECPurdkAOD+vcZTxq+N1mPUicsBPAHng0NAffuZ+poudR1HWUN9daC+pktdx1EeWKWzI6NIbgnenHzis42Y060//fq2WxkJ2F5pROQSOI+IjlbVQc/uuvYz9TVd6jqOUof66gv1NV3qOo4ywUadLZMnei2cOBjvHcEB6Lx7IN1Z5f4O689VAA5ofTvV/Xu0p4xfG0CNrouIXArg8wA+oqqLWnbVtZ+pr+lS13GUCdTXDqiv6VLXcZQZtupsaYxoVe0DMA/AOM+ucWiPcyFmLIYzaJr9KSK7AzgGw/05G04Q/wdb6n0QwF6eMse4dRuMg/Nm7JIsBLcNEbkcwKlwlPs5z+5a9jP1NXVqOY6ygPraCfU1dWo5jrLCap0t+k3LiG9lngKgD8DpcN62vBxOsPghRctm4487aN7j/mwH8B3374Pd/f8N523Uz8BJ6fIb+KeFmQ/gA+6gm4/2tDCvdgfwb9w2PuO2WYsUPAAmuuf7ETh3sY2fvVvK1LKfqa+R+4v6mn0fU1+D+4b6Gq2/qK/59LPVOlt4B8Xo0LPg3BXshHPn/OGiZbL1B8BxcGJ5vD+T3f0C4AI46WF6ATwE4J2eNvYDcLM7mDa7f+/rKXMEnDdhe922zkdN0u8E9K8CuKClTG37mfoaqa+or9n3MfU1vH+or+Z9RX3Np5+t1llxKxJCCCGEEEIMKU1MNCGEEEIIIbZAI5oQQgghhJCI0IgmhBBCCCEkIjSiCSGEEEIIiQiNaEIIIYQQQiJCI5oQQgghhJCI0IgmhBBCCCEkIjSiCSGEEEIIiQiNaEIIIYQQQiLy/wE5lem7yYYSKQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x216 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plotXs(x_true, x_rml, x_rnd, 'fault.pdf')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
